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Abstract 

Spin-orbit coupling can be described in two approaches. The first method, known as the 
"MacDonald torque", is often combined with a convenient assumption that the quality factor 
Q is frequency-independent. This makes the method inconsistent, for the MacDonald theory 
tacitly fixes the rheology of the mantle by making Q scale as the inverse tidal frequency. 

Spin-orbit coupling can be treated also in an approach called "the Darwin torque". 
While this theory is general enough to accommodate an arbitrary frequency-dependence 
of Q , this advantage has not yet been fully exploited in the literature, where Q is often 
assumed constant or is set to scale as inverse tidal frequency, the latter assertion making 
the Darwin torque equivalent to a corrected version of the MacDonald torque. 

However neither a constant nor an inverse-frequency Q reflect the properties of realistic 
mantles and crusts, because the actual frequency-dependence is more complex. Hence it is 
necessary to enrich the theory of spin-orbit interaction with the right frequency-dependence. 

We accomplish this programme for the Darwin-torque-based model near resonances. We 
derive the frequency-dependence of the tidal torque from the first principles of solid-state 
mechanics, i.e., from the expression for the mantle's compliance in the time domain. We 
also explain that the tidal torque includes not only the customary, secular part, but also an 
oscillating part. 

We demonstrate that the Impq term of the Darwin-Kaula expansion for the tidal torque 
smoothly passes zero, when the secondary traverses the Impq resonance (e.g., the principal 
tidal torque smoothly goes through nil as the secondary crosses the synchronous orbit). 

Thus we prepare a foundation for modeling entrapment of a despinning primary into a 
resonance with its secondary. The roles of the primary and secondary may be played, e.g., 
by Mercury and the Sun, correspondingly, or by an icy moon and a Jovian planet. 

We also offer a possible explanation for the unexpected frequency-dependence of the 
tidal dissipation rate in the Moon, discovered by LLR. 

1 Introduction 

We continue a critical examination of the tidal-torque techniques, begun in Efroimsky & Williams 
(2009), where the empirical treatment by and MacDonald (1964) was considered from the view- 
point of a more general and rigorous approach by Darwin (1879, 1880) and Kaula (1964). Re- 
ferring the Reader to Efroimsky & Williams (2009) for proofs and comments, we begin with an 
inventory of the key formulae describing the spin-orbit interaction. While in Ibid, we employed 
those formulae to explore tidal despinning well outside the 1:1 resonance (and in neglect of the 
intermediate resonances) , in the current paper we apply this machinery to the case of despinning 
in the vicinity of a spin-orbit resonance. 
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Although the topic has been popular since mid-sixties and has already been addressed in 
books, the common models are not entirely adequate to the actual physics. Just as in the 
nonresonant case discussed in Ibid., a generic problem with the popular models of libration or of 
capture into a resonance is that they employ wrong rheologies (the work by Rambaux et al. 2010 
being the only exception we know of). Above that, the model based on the MacDonald torque 
suffers a defect stemming from a genuine inconsistency inherent in the theory by MacDonald 
(1964). 

As explained in Efroimsky and Williams (2009) and WiUiams and Efroimsky (2012), the 

MacDonald theory, both in its original and corrected versions, tacitly fixes an unphysical shape 
of the functional dependence Q{x) j where Q is the dissipation quality factor and x is the tidal 
frequency (Williams & Efroimsky 2012). So we base our approach on the developments by Darwin 
(1879, 1880) and Kaula (1964), combining those with a reahstic law of frequency- dependence of 
the damping rate. 

Since our main purpose is to lay the groundwork for the subsequent study of the process of 
falling into a resonance, the two principal results obtained in this paper are the following: 

(a) Starting with the realistic rheological model (the expression for the compliance in the time 
domain) , we derive the complex Love numbers as functions of the frequency x > write down 
their negative imaginary parts as functions of the frequency: —Im[ki{x)] — |^/(x)| sine^(x). 
It is these expressions that appear as factors in the terms of the Darwin-Kaula expansion of tides. 
These factors' frequency-dependencies demonstrate a nontrivial shape, especially near resonances. 
This shape plays a crucial role in modeling of despinning in general, specifically in modeling the 
process of falling into a spin-orbit resonance. 

(b) We demonstrate that, beside the customary secular part, the Darwin torque contains a 
usually omitted oscillating part. 

2 Linear bodily tides 

Linearity of tide means that: (a) under a static load, deformation scales linearly, and (b) under 
undulatory loading, the same linear law applies, separately, to each frequency mode. The latter 
implies that the deformation magnitude at a certain frequency should depend linearly upon the 
tidal stress at this frequency, and should bear no dependence upon loading at other tidal modes. 
Thence the dissipation rate at that frequency will depend on the stress at that frequency only. 

2.1 Linearity of the tidal deformation 

— ♦ 

At a point R = (i?. A, 0), the potential due to a tide-raising secondary of mass M*^^, located at 
r* — (r*. A*, (p*) with r* > R,is expandable over the Legendre polynomials P;(cos7) : 

I 

1=2 1=2 ^ ^ 

oo ^ I 

^ 1=2 ' m=0 ^ 

where G = 6.7 x lO^^"*^ m'^kg^^s^^ is Newton's gravity constant, and 7 is the angular separation 
between the vectors f* and R pointing from the primary's centre. The latitudes 0, 0* are reck- 
oned from the primary's equator, while the longitudes A, A* are reckoned from a fixed meridian. 

Under the assumption of linearity, the Z term Wi{R, r*) in the secondary's potential 
causes a linear deformation of the primary's shape. The subsequent adjustment of the primary's 
potential being linear in the said deformation, the / adjustment Ui of the primary's potential 
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is proportional to Wi . The theory of potential requires Ui{r) to fall off, outside the primary, as 
^-{i+i) 'jj^ug the overall amendment to the potential of the primary amounts to: 

1=2 1=2 ^ ^ 

R now being the mean equatorial radius of the primary, R = {R, (p, X) being a surface point, 
r = (r , , A) being an exterior point located above it at a radius r > R. The coefficients ki , 
called Love numbers, are defined by the primary's rheology. 

For a homogeneous incompressible spherical primary of density p , surface gravity g, and 
rigidity fi , the static Love number of degree / is given by 

3 1 , ^ (2/2 + 4/ + 3) 3 (2/2 + 4/ + 3)// 
= W^) TTT, ' "^^^^^ " J^;^^ = AlrrGp^R^ " 

For R <^ r , r* , consideration of the I = 2 input in ([2]) turns out to be sufficient 

These formulae apply to static deformations. However an actual tide is never static, except 
in the case of synchronous orbiting with a zero eccentricity and inclinationj^l Hence a realistic 
perturbing potential produced by the secondary carries a spectrum of modes w^^pg (positive 
or negative) numbered with four integers Impq as in formula fllOSp below. The perturbation 
causes a spectrum of stresses in the primary, at frequencies Ximpq ~ \^impq \ ■ Although in a 
linear medium strains are generated exactly at the frequencies of the stresses, friction makes 
each Fourier component of the strain fall behind the corresponding component of the stress. 
Friction also reduces the magnitude of the shape response - hence the deviation of a dynamical 
Love number ki{x) from its static counterpart ki = ki{0) . Below we shall explain that formulae 
([2] -[3]) can be easily adjusted to the case of undulatory tidal loads in a homogeneous planet or 
in tidally-despinning homogeneous satellite (treated now as the primary, with its planet playing 
the role of the tide-raising secondary) . However generalisation of formulae ([2] - [3]) to the case of 
a librating moon (treated as a primary) turns out to be highly nontrivial. As we shall see, the 
standard derivation by Love (1909, 1911) falls apart in the presence of the non-potential inertial 
force containing the time- derivative of the primary's angular velocity. 

The frequency-dependence of a dynamical Love numbers takes its origins in the "inertia" of 
strain and, therefore, of the shape of the body. Hence the analogy to linear circuits: the /*^ 
components of W and U act as a current and voltage, while the /*^ Love number plays, up to a 
factor, the role of impedance. Therefore, under a sinusoidal load of frequency x 5 it is convenient 
to replace the actual Love number with its complex counterpart 

^lix) = I ^/(x) I exp [- ze/(x)] , (4) 

ei being the frequency- dependent phase delay of the reaction relative to the load (Munk & 
MacDonald 1960, Zschau 1978). The "minus" sign in (j4]) makes U lag behind W for a positive 
ei . (So the situation resembles a circuit with a capacitor, where the current leads voltage.) 

In the limit of zero frequency, i.e., for a steady deformation, the lag should vanish, and so 
should the entire imaginary part: 

Im [kiO)] = \ki{0)\ sinQ(O) = , (5) 



^ Special is the case of Phobos, for whose orbital evolution the and perhaps even the ^4 terms may be 
relevant (Bills et al. 2005). Another class of exceptions is constituted by close binary asteroids. The topic is 
addressed by Taylor & Margot (2010), who took into account the Love numbers up to fcg . 

^ The case of a permanently deformed moon in a 1:1 spin-orbit resonance falls under this description too. 
Recall that in the tidal context the distorted body is taken to be the primary. So from the viewpoint of the 
satellite its host planet is orbiting the satellite synchronously, thus creating a static tide. 
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leaving the complex Love number real: 



kiO) = 7^e[^^(0)] = I fc; (0)1 cos 6^(0) , (6) 

and equal to the customary static Love number: 

ki{0) = ki . (7) 

Solution of the equation of motion combined with the constitutive (rheological) equation 
renders the complex ki{x)i explained in Appendix lD.il Once ki{x) is found, its absolute value 

kiix) = \ki{x)\ (8) 

and negative argument 

[ki{x)\ 

should be inserted into the term of the Fourier expansion for the tidal potential. Things get 
simplified when we study how the tide, caused on the primary by a secondary, is acting on that 
same secondary. In this case, the l^^ term in the Fourier expansion contains ^^nd ^i{x) 

in the convenient combination ki{x) sine;(x) , which is exactly — Xm . 

Rigorously speaking, we should say not "the l^^ term", but "the l^^ terms", as each / cor- 
responds to an infinite set of positive and negative Fourier modes uimpq , the physical forcing 
frequencies being x = Ximpq = \<^impq\ ■ Thus, while the functional forms of both and 
sine/(x) depend only on / , both functions take values that are different for different sets of num- 
bers mpq. This happens because x assumes different values Ximpq on these sets. Mind though 
that for triaxial bodies the functional forms of and sine;(x) may depend also on m, p, q. 

2.2 Damping of a linear tide 

Beside the standard assumption Ui{r) oc Wi{R, r*) , the linearity condition includes the re- 
quirement that the functions ki{x) and e^(x) be well defined. This implies that they depend 
solely upon the frequency Xi and not upon the other frequencies involved. Nor shall the Love 
numbers or lags be influenced by the stress or strain magnitudes at this or other frequencies. 

Then, at frequency the mean (over a period) damping rate {E (x)) depends on the value 
of X and on the loading at that frequency, and is not influenced by the other frequencies: 

{E{x)) = -^^^^ (10) 

Qix) 

or, equivalently: 

AE^ycUx) = - ^^^1^ ^ (11) 

AEcydeix) being the one-cycle energy loss, and Q{x) being the so-called quality factor. 

If Epeakix) ill dini - [HD is agreed to denote the peak energy stored at frequency Xi the 
appropriate Q factor is connected to the phase lag e(x) through 

Q-^ = sin lei . (12) 

^ energy ' ' ^ ■' 

and not through Q^^^^^y =tan|e| as often presumed (see Appendix [B] for explanation). 
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If Epeakix) is defined as the peak work, the corresponding Q factor is related to the lag via 



tan 


e 






e 


) 


tan 


e 



^ work / „ . . \ . . ' (■^'^) 



as demonstrated in Appendix [B] below El In the limit of a small e , (fT3|) becomes 

Ql\^^ = sin|6| + 0{e') = \e\ + 0{e') , (14) 

so definition (fT3|) makes 1/Q a good approximation to sin e for small lags only. 

For the lag approaching 7r/2, the quality factor defined through (IT^ attains its minimum, 
Q energy ~ ^ ' whilc definition (IT^ furnishes Q^^i. = . The latter is not surprising, as in the 
said limit no work is carried out on the system. 

Linearity requires the functions kiix) ^'^d therefore also ei{x) to be well-defined, i.e., to be 
independent from all the other frequencies but x- We now see, the requirement extends to Q{x) ■ 

The third definition of the quality factor (offered by Golderich 1963) is Q^^^i^^^^^j.^ = tan |e| . 
However this definition corresponds neither to the peak work nor to the peak energy. The 
existing ambiguity in definition of Q makes this factor redundant, and we mention it here only as 
a tribute to the tradition. As we shall see, all practical calculations contain the products of the 
Love numbers by the sines of the phase lags, ki sin ei , where / is the degree of the appropriate 
spherical harmonic. A possible compromise between this mathematical fact and the historical 
tradition of using Q would be to define the quality factor through (IT^ . in which case the quality 
factor must be equipped with the subscript / . (This would reflect the profound difference between 
the tidal quality factors and the seismic quality factor - see Efroimsky 2012.) 

3 Several basic facts from continuum mechanics 

This section offers a squeezed synopsis of the basic facts from the linear solid-state mechanics. 
A more detailed introduction, including a glossary and examples, is offered in Appendix |Al 

3.1 Stationary linear deformation of isotropic incompressible media 

Mechanical properties of a medium are furnished by the so-called constitutive equation or con- 
stitutive law, which interrelates the stress tensor S with the strain tensor U defined as 



2 



(V ® u) + (V ® u) , (15) 



where u is the vector of displacement. 

As we shall consider only linear deformations, our constitutive laws will be linear, and will 
be expressed by equations which may be algebraic, differential, integral, or integro-differential. 

The elastic stress § is related to U through the simplest constitutive equation 

§ = B U , (16) 

B being a four-dimensional matrix of real numbers called elasticity moduli. 

{h) 

A hereditary stress § is connected to U as 

(h) 

§ = B U , (17) 



■^Deriving this formula in Appendix to Efroimsky & Williams (2009), we inaccurately termed £'peafc(x) 
peak energy. However our calculation of Q was carried out in understanding that Epeak (x) is the peak work. 
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B being a four- dimensional integral-operator-valued matrix. Its component Bij^i acts on an 
element um of the strain not as a mere multiplier but as an integral operator, with integration 
going from t' = — oo through t' = t . To furnish the value of CTjj = Yliki BijkiUki at time t , the 
operator "consumes" as arguments all the values of Uki{t') over the interval t' G ( — oo, t ] . 
The viscous stress is related to the strain through a differential operator A ^ : 

{^) d 

§ = A - U , (18) 
ot 

A being a four-dimensional matrix consisting of empirical constants called viscosities. 

In an isotropic medium, each of the three matrices, B , B , and A , includes two terms only. 
The elastic stress becomes: 



S = sLnrnet^c + ,e.^ator^c = 3KQlSpU^+2/^(^U-ilSpU^ , (19) 



with K and fi being the bulk elastic modulus and the shear elastic modulus, correspondingly, I 
standing for the unity matrix, and Sp denoting the trace of a matrix: Sp U = J2i ■ 
The hereditary stress becomes: 



{h) {h) (h) ~ f I \ f I \ 

S = ^.olumetr^c+ § = 3 - 1 Sp U j + 2 /i U - - 1 Sp U j , (20) 



where K and are the hulk-modulus operator and the shear-modulus operator, accordingly. 
The viscous stress acquires the form: 

b') W (f) d fl \ d f I \ 

S = ^.olumetr^c + § de«a*o„c = ^ C ^ (^3 I Sp U j + 2 7/ - U - - I Sp U j , (21) 

the quantities C, and rj being termed as the hulk viscosity and the shear viscosity, correspondingly 

1 1 
The term ^ I Sp U is called the volumetric part of the strain, while U — ^ I Sp U is called 

the deviatoric part. Accordingly, in expressions (11611 - [T63ll for the stresses, the pure-trace terms 

are called volumetric, the other term being named deviatoric. 

If an isotropic medium is also incompressible, the relative change of the volume vanishes: 

Sp U = , and so does the expansion rate: V-v = ^SpU = 0. Then the volumetric part of 

the strain becomes zero, and so do the volumetric parts of the elastic, hereditary, and viscous 

stresses. The incompressibility assumption may be applicable both to crusty objects and to large 

icy moons of low porosity. At least for lapetus, the low-porosity assumption is likely to be correct 

(Castillo-Rogez et al. 2011). 

3.2 Approaches to modeling viscoelastic deformations. 
Problems with terminology 

One approach to linear deformations is to assume that the elastic, hereditary and viscous devia- 
toric stresses simply sum up, each of them being linked to the same overall deviatoric strain: 

{total) (e) {h) [v) , d f ^ d \ 

S = S + S + S = BU+BU + A— U= B + B + A— U . 22 

ot \ ot J 

An alternative option, to be used in section 15.31 below, is to start with an overall deviatoric 
stress, and to expand the deviatoric strain into elastic, viscous, and hereditary parts: 

(e) {h) {v) (e) 1 {v) 1 /•* (h) 

U=U + U + U , U=-S, U=- / ^{t')dt' , u= , (23) 
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J being an integral operator with a time-dependent kernel. 

An even more general option would be to assume that both the strain and stress are comprised 
by components of different nature - elastic, hereditary, viscous, or more comphcated (plastic). 
Which option to choose - depends upon the medium studied. The rich variety of materials offered 
to us by nature leaves one no chance to develop a unified theory of deformation. 

As different segments of the continuum-mechanics community use different conventions on 
the meaning of some terms, we offer a glossary of terms in Appendix El Here we would only 
mention that in our paper the term viscoelastic will be applied to a model containing not only 
viscous and elastic terms, but also an extra term responsible for an anelastic hereditary reaction. 
(A more appropriate term viscoelastohereditary would be way too cumbersome.) 

3.3 Evolving stresses and strains. Basic notations 

In the general case, loading varies in time, so one has to deal with the stress and strain tensors as 
functions of time. However, treatment of viscoelasticity turns out to be simpler in the frequency 
domain, i.e., in the language of complex rigidity and complex compliance. To this end, the stress 
cr-yy and strain u^^ in a linear medium can be Fourier-expanded as 



n=0 



n=0 



iXj^t + ifa- (Xn) 



(24a) 



E 

n=0 



Tie 



(24b) 



(t) = J2 U^uiXn) COs[xJ + MXn)] = J] 



n=0 



n=0 



U^AXr. 



(25a) 



J2 u^uiXn) e 



n=0 



(25b) 



where the complex amplitudes are: 



(26) 



while the initial phases fa{x) V^u(x) ^^^e chosen in a manner that sets the real amplitudes 
(^-lAXn) and u^uiXn) non-negative. 

We wrote the above expansions as sums over a discrete spectrum, as the spectrum generated 
by tides is discrete. Generally, the sums can, of course, be replaced with integrals over frequency: 



ixt 

(^•yvix) 6 dx and u^u{t) 



ixt 



u-yuix) e dx 



(27) 



Whenever necessary, the frequency is set to approach the real axis from below: Xm{x) — )■ 0- 
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3.4 Should we consider positive frequencies only? 

At first glance, the above question appears pointless, as a negative frequency is a mere abstraction, 
while physical processes go at positive frequencies. Mathematically, a full Fourier decomposition 
of a real field can always be reduced to a decomposition over positive frequencies only. 
For example, the full Fourier integral for the stress can be written as 



iuit 



ixt — ix* 



(28) 



where we define x = \^ \ ■ Denoting complex conjugation with asterisk, we write: 



Kui-X) e + s*{x) e 



dx . 



(29) 



The stress is real: o"* (t) = a^y{t) . Equating the right-hand sides of (|28|1 and (|29|) . we obtain 



"^71^ ( X) ^-yu^X) 



(30) 



whence 



ixt — ixt 

s^uix) e + s*{x) e 



/ _ ixt 

dx = 2 s^^(x) e dx ■ (31) 

Jo 



This leads us to fl271). if we set 



a. 



ix) = '^Sjuix) 



(32) 



While the switch from 



iuit 



"yuit) = s^u^Lo) e duj to the expansion a^y{t) = a^^i 
makes things simpler, the simplification comes at a cost, as we shall see in a second. 
Recall that the tide can be expanded over the modes 



CO 



ixt 

e dx 



ujimpq = {I -2p) CO + {I -2p + q) M + m {tl 



2p + q)n — mO 



(33) 



each of which can assume positive or negative values, or be zero. Here /, m, q are some integers, 
9 is the primary's sidereal angle, 9 is its spin rate, while w, fi, Ai and n are the secondary's 
periapse, node, mean anomaly, and mean motion. The appropriate tidal frequencies, at which 
the medium gets loaded, are given by the absolute values of the tidal modes: Xi 



mpq I ^Impq 



The positively-defined forcing frequencies Ximpq actual physical frequency at which 

the Impq term in the expansion for the tidal potential (or stress or strain) oscillates. 

The motivation for keeping also the modes w^^p^ is subtle: it depends upon the sign of w^.^^^ 
whether the Impq component of the tide lags or advances. Specifically, the phase lag between 
the Impq component of the perturbed primary's potential U and the Impq component of the 
tide-raising potential W generated by the secondary is given by 



'Impq 



^Impq ^^Impq 



Impq 



At 



Impq 



sgncu 



Impq 



Y At 

^ Impq Impq 



sgnu 



Impq 



(34) 



where the time lag Ati^^^ is always positive. 

While the lag between the applied stress and resulting strain in a sample of a medium is 
always positive, the case of tides is more complex: there, the lag can be either positive or 
negative. This, of course, in no way implies whatever violation of causality (the time lag At^^^^ 
is always positive). Rather, this is about the directional difference between the planetocentric 
positions of the tide-raising body and the resulting bulge. For example, the principal component 
of the tide, Impq = 2200 , stays behind (has a positive phase lag 62200 ) when the secondary is 



below the synchronous orbit, and advances (has a negative phase lag 



^2200 



when the secondary 
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is at a higher orbit. To summarise, decomposition of a tide over both positive and negative 
modes (^i^^pq {8Jid not just over the positive frequencies Ximpq ) ^^^^ have a physical meaning, 
as the sign of a mode carries physical information. 

Thus we arrive at the following conclusions: 



1. As the fields emerging in the tidal theory - the tidal potential, stress, and strain - are all 
real, their expansions in the frequency domain may, in principle, be written down using the 
positive frequencies x only. 



In the tidal theory, the potential (and, consequently, the tidal torque and force) contain 
components corresponding to the tidal modes w^^p^ of both the positive and negative signs. 
While the Impq components of the potential, stress, and strain oscillate at the positive 
frequencies Ximpq ~ \^impq\ ' ^^S^ each uj^^^^ does carry physical information: it dis- 
tinguishes whether the lagging of the Impq component of the bulge is positive or negative 
(falling behind or advancing). Accordingly, this sign enters explicitly the expression for the 
appropriate component of the torque or force. Hence a consistent tidal theory should be 
developed through expansions over both positive and negative tidal modes cj^^^g and not 
just over the positive Ximpq ■ 

In order to rewrite the tidal theory in terms of the positively-defined frequencies Ximpq 
only, one must inserts "by hand" the extra multipliers 



{l-2p + q)n - mO (35) 



into the expressions for the Impq components of the tidal torque and force. 



4. One can employ a rheological law (constitutive equation interconnecting the strain and 
stress) and a Navier- Stokes equation (the second law of Newton for an element of a vis- 
coelastic medium), to calculate the phase lag e^^^^ of the primary's potential Ui^^^ relative 
to the potential W^^^^ generated by the secondary. If both these equations are expanded, 
in the frequency domain, via positively-defined forcing frequencies Ximpq '^'^^y^ the resulting 
phase lag, too, will emerge as a function of Ximpq '■ 

^Impq = ^liXlmpq) ■ (36) 

Within this treatment, one has to equip the lag, "by hand", with the multiplier (135|) . 



As we saw above, the lag (l36l) is the argument of the complex Love number k^^Ximpq) ■ Solution 
of the constitutive and Navier-Stokes equations renders the complex Love numbers, from which 
one can calculate the lags. Hence the above item [4] may be rephrased in the following manner: 

4'. Under the convention that f/^^^^ = ?7(x;„pJ and 1^^^^^ = IV(Xi^pg) , we have: 

Ulmpq = kiXlmpq)^lmpq ^heU W^^p, > , i.e., whcU W^^^^ = Ximpq ' (37a) 

Ulmpq = K(Xlmpq)^lmpq ^hcU W^^p^ < , i.e., whcU UJ impq Ximpq ^ (37b) 

asterisk denoting the complex conjugation. 
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This ugly convention, a switch from to k* , is the price we pay for employing only the positive 
frequencies in our expansions, when solving the constitutive and Navier-Stokes equations, to find 
the Love number. In other words, this is a price for our pretending that Wi^^^ and f/^^p^ are 
functions of Ximpq ^ whereas in reahty they are functions of w^^pg • 

Alternative to this would be expanding the stress, strain, and the potentials over the positive 
and negative modes cu^^pg , with the negative frequencies showing up in the equations. With the 
convention that f/^^p^ = t/(w,^pg) and W^^p^ = Vr(w,^p,) , we would have 

Ulmpq = ki^lmpq)^lmpq ' ^^^11 W^^p^ . (38) 



All these details can be omitted at the despinning stage, if one keeps only the leading term of 
the torque and ignores the other terms. Things change, though, when one takes these other terms 
into account. On crossing of an Impq resonance, factor (1351) will change its sign. Accordingly, 
the Impq term of the tidal torque (and of the tidal force) will change its sign too. 



3.5 The complex rigidity and compliance. Stress-strain relaxation 

The stress cannot be obtained by means of an integral operator that would map the past history 
of the strain, U(t') over t' G ( — oo, t ] , to the value of § at time t. The insufficiency of 
such an operator is evident from the presence of a time-derivative on the right-hand side of (ITSl) . 
Exceptional are the cases of no viscosity (e.g., a purely elastic material). 

On the other hand, we expect, on physical grounds, that the operator J inverse to /i is an 
integral operator. In other words, we assume that the current value of the strain depends only 
on the present and past values taken by the stress and not on the current rate of change of the 
stress. This assumption works for weak deformations, i.e., insofar as no plasticity shows up. So 
we assume that the operator J mapping the stress to the strain is just an integral operator. 

Since the forced medium "remembers" the history of loading, the strain at time t must 
be a sum of small installments ^ ^(^ ~ ^') da^iy{t') , each of which stems from a small change 
da^u{t — r) of the stress at an earlier time t' < t. The entire history of the past loading results, 
at the time t , in a total strain Uyt,{t) rendered by an integral operator J{t) acting on the entire 
function a^^it') and not on its particular value (Karato 2008): 

2uj^{t) = J{t)a^u = / Jir) a^^{t-T)dT = J{t - t') a^^ {t') dt' , (39) 

where t ' is some earlier time ( t ' < t), overdot denotes d/dt' , while the "age variable" t = t — t' 
is reckoned from the current moment t and is aimed back into the past. The so-defined integral 
operator J{t) is called the compliance operator, while its kernel J{t — t') goes under the name 
of the compliance function or the creep-response function. 

Integrating fl39l) by parts, we recast the compliance operator into the form of 

POO 

2u^y{t) = J{t) a^y = J{Q) a^yif) — J(oo) cr^^(— oo) + / J {t) a^^it — r) dr (40a) 



= J{0)ajy{t) - J(oo) (T^^(-oo) + f j {t - t') a^y{t') dt' . (40b) 

J — oo 

The quantity J(oo) is the relaxed compliance. Being the asymptotic value of J{t — t') at 
t — t'^ oo , this parameter corresponds to the strain after complete relaxation. The load in the 
infinite past may be assumed zero, and the term — J(oo) cr^i.(— oo) may be dropped 
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The second important quantity emerging in ( HOj) is the unrelaxed compliance J(0) , which is 
the value of the comphance function J{t — t') at t — t' = 0. This parameter describes the 
instantaneous reaction to stressing, and thus defines the elastic part of the deformation (the 
rest of the deformation being viscous and hereditary). Thus the term containing the unrelaxed 
compliance J(0) should be kept. The term, though, can be absorbed into the integral if we agree 
that the elastic contribution enters the compliance function not as 

J{t — t') = J(0) + viscous and hereditary terms , (41) 

but as 

J(t — t') = J(0) 0(t — t') + viscous and hereditary terms , (42) 

the Heaviside step- function 6(t — t') being unity for t — t' > , and zero for t — t' < . As 
the derivative of the step-function is the delta-function 6{t — t ') , we can write ( ]40bl) simply as 

2u^y{t) = J{t) a^y = I j {t-t') a^^{t')dt' , with J{t-t') containing J(O)0(t-t') 

J — oo 

Equations ( 139|) . ( l40l) . ( l43l) are but different expressions for the compliance operator J acting as 

2u^y = J ■ (44) 

Inverse to the compliance operator is the rigidity operator jl defined through 

a^u = 2 /t u^fu . (45) 

Generally, /t is not just an integral operator, but is an integro-differential operator. So it cannot 
take the form of (T^y{t) = 2 /i(t — t') u^^{t') dt' . However it can be written as 

a^^(t) = 2 I fi{t - t')u^^{t')dt' , (46) 

J — oo 

if we permit the kernel fi{t — t') to contain a term r/ 5(t—t '), where 6{t—t') is the delta-function. 
After integration, this term will furnish the viscous part of the stress, 2 rj ii^^ . 

The kernel ^{t — t') goes under the name of the stress-relaxation function. Its time- 
independent part is n{0) Q{t — t') , where the unrelaxed rigidity /i(0) is inverse to the unrelaxed 
compliance J(0) and describes the elastic part of deformation. Each term in /i(t — t') , which 
neither is a constant nor contains a delta-function, is responsible for hereditary reaction. 

For more details on the stress-strain relaxation formalism see the book by Karato (2008). 

3.6 Stress-strain relaxation in the frequency domain 

Let us introduce the complex compliance J{x) and the complex rigidity fi{x) , which are, by 
definition, the Fourier images not of the J(r) and /u(t) functions, but of their time-derivatives|j 

/•oo /»oo 

/ J(x)e»^-dx = j(r) , where J(x) = / j (r) e'^^^ dr . (47) 
Jo Jo 



^ Expressing the stress through the strain, we encountered three possibihties: the elastic stress was simply 
proportional to the strain, the viscous stress was proportional to the time-derivative of the strain, while the 
hereditary stress was expressed by an integral operator jj . However, when we express the strain through the 
stress, we place the viscosity into the integral operator, so the purely viscous reaction also looks like hereditary. 
It is our convention, though, to apply the term hereditary to delayed reactions other than purely viscous. 

^ Recall that it is the time-derivative of J(t) that is the kernel of the integral operator (j43|) . Hence, to arrive 

at ([50l) , we have to define J(x) as the Fourier image of J (r) . 
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and 



POO POO 

/ fi{x)e'''^dx = Kt) , where fiix) = / A(r)e-^^"rfr , (48) 
Jo Jo 

the integrations over r spanning the interval [ 0, oo ) , as both kernels are nil for r < anyway. 
In P7|) and (HHj) . we made use of the fact (explained in subsection I3.4p that, when expanding 
real fields, it is sufficient to use only positive frequencies. 

Expression (139|1 . in combination with the Fourier expansions (127|) and with (H7|) . furnishes: 

roo /-oo 

2 / y-yuix) e dx = ^iiuix) J{x) e rfx , (49) 
Jo Jo 

which leads us to: 

"^Uyuix) = Jix)^yAx) ■ (50) 

Similarly, insertion of f l27|) into f H6|) leads to the relation 

0-7:.(x) = 2 /i(x) M7;.(X) , (51) 

comparison whereof with ( !50|) immediately entails: 

J(X) /i(x) = 1 . (52) 
Writing down the complex rigidity and compliance as 

fiix) = \fi{x)\ exp[z5(x)] (53) 

and 

Ax) = \Ax)\ exp[-z5(x)] , (54) 
we split (|52|) into two expressions: 

I^(X)I = (55) 

Hx)\ 

and 



¥>u{x) = ¥><t{x) - Six) ■ (56) 

From the latter, we see that the angle 6{x) is a measure of lagging of a strain harmonic mode 
relative to the appropriate harmonic mode of the stress. It is evident from (13^ - that 

tan^fy) = - ^"^^^^^^^ - ^^^'^^^^^ (57) 



4 Complex Love numbers 

The developments presented in this section will rest on a very important theorem from solid- 
state mechanics. The theorem, known as the correspondence principle, also goes under the name 
of elastic-vis coelastic analogy. The theorem applies to linear deformations in the absence of 
nonconservative (inertial) forces. While the literature attributes the authorship of the theorem 
to different scholars, its true pioneer was Sir George Darwin (1879). One of the corollaries ensuing 
from this theorem is that, in the frequency domain, the complex Love numbers are expressed via 
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the complex rigidity or compliance in the same way as the static Love numbers are expressed via 
the relaxed rigidity or compliance. 

As was pointed out much later by Biot (1954, 1958), the theorem is inapplicable to non- 
potential forces. Hence the said corollary fails in the case of librating bodies, because of the 
presence of the inertial forc^ —C3 x rp, where p is the density and uj is the libration angular 
velocity. So the standard expression (j3]) for the Love numbers, generally, cannot be employed for 
librating bodies. 

Subsection 14.11 below explains the transition from the stationary Love numbers to their dy- 
namical counterparts, the so-called Love operators. We present this formalism in the frequency 
domain, in the spirit of Zahn (1966) who pioneered this approach in application to a purely vis- 
cous medium. Subsection 14.21 addresses the negative tidal modes emerging in the Darwin-Kaula 
expansion for tides. Employing the correspondence principle, in subsection 14.31 we then write 
down the expressions for the factors |fc;(x)| sine^(x) = "^^"^[^/(x)] emerging in the expansion 
for tides. Some technical details of this derivation are discussed in subsections 14.41 - H75l 

For more on the correspondence principle and its applicability to Phobos see Appendix 

4.1 From the Love numbers to the Love operators 

A homogeneous incompressible primary, when perturbed by a static secondary, yields its form 
and, consequently, has its potential changed. The l^^ spherical harmonic Ui{r) of the resulting 
increment of the primary's exterior potential is related to the l*^ spherical harmonic Wi{R,r) 
of the perturbing exterior potential through ([2]). 

As the realistic disturbances are never static (except for synchronous orbiting), the Love 
numbers become operators: 



Ui{r,t)=[^-j ki{t)Wi{R, r*,t') . (58) 

A Love operator acts neither on the value of W at the current time t , nor at its value at an 
earlier time t' , but acts on the entire shape of the function Wi{R, r*, t') , with t' belonging 
to the semi-interval (— c>o, t) . This is why we prefer to write ki{t) and not ki{t, t') . 
Being linear for weak forcing, the operators must read: 

- / ki{T)WiiR,r*,t-r)dr= (-] / ki{t - t') Wi {R , r* , t') dt' (59^) 

or, after integration by parts: 

-j [k{0)W{t) - ki{oo)W{-oo)]+ (-] kir) Wi{R,r\t-r)dT (59b) 



r 



[ki{Q)W{t) - ki{oo)W{-oo)] + j k,{t-t') Wi{R,r*,t')dt' (59c) 



- ki(oo)W(-oo) 
r J 

7) ^ / ^'^^ " ^'^ " + " ^'^ ] Wi{R,r*,t')dt' . (59d) 



^ The centripetal term is potential and causes no troubles, except for the necessity to introduce a degree-0 
Love number. 
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Just as in the case of the comphance operator (1391- 1401) . in expressions ( l59l) we obtain the terms 
ki{0)W(t) and — ki{oo)W{—oo) . Of the latter term, we can get rid by setting W{—oo) nil, 
while the former term may be incorporated into the kernel in exactly the same way as in (jUj- 
l^3|) . Thus, dropping the unphysical term with W{—oo) , and inserting the elastic term into the 
Love number not as ki{0) but as ki{0) Q(t — t') , we simplify (159d|) to 



Ui{r,t) = (^^Y j' k,{t-t')Wi{R,r*,t')dt' , (60) 

with ki{t — t') now including, as its part, ki{0)Q(t — t') instead of ki{0) . 

Were the body perfectly elastic, ki{t — t') would consist of the instantaneous- reaction term 
ki{0)Q{t — t') only. Accordingly, the time-derivative of ki would be: k^lt — t') = ki6{t — t') 
where ki = ki{0) , so expressions (jSH-En]) would coincide with ([2]). 

Similarly to introducing the complex compliance, one can define the complex Love numbers 

as Fourier transforms of ki (t) : 

r>00 

kix)e''''dx =k (r) , (61) 







the overdot standing for d/dr . Churkin (1998) suggested to term the time-derivatives ki {t) 
as the Love functions^ Inversion of f l6T]) trivially yields: 



kiix) = / Mr) e-'^^dr = k,{0) + zx / ^r) - A:,(0)e(r)] e^^^^ rfr , (62) 
Jo Jo 

where we integrated only from because the future disturbance contributes nothing to the 
present distortion, so ki{T) vanishes at r < 0. Recall that the time r denotes the difference 
t — t'. So r is reckoned from the present moment t and is directed back into the past. 

Defining in the standard manner the Fourier components Uiix) and W^lx) of functions 
Ui (t) and Wi (t) , we write (l59l) in the frequency domain: 



Uiix) = k,{x) W,{x) , (63) 

where we denote the frequency simply by x instead of the awkward Ximpq ■ employ ( 163|) in the 
tidal theory, one has to know the frequency- dependencies k^ (x) . 



4.2 The positive forcing frequencies x = l*^! vs. 
the positive and negative tidal modes u 

It should be remembered that, by relying on formula (!63|) . we place ourselves on thin ice, because 
the similarity of this formula to ( 150|) and ( 15T|) is deceptive. 

In ( 15U]1 and (ICTl . it was legitimate to limit our expansions of the stress and the strain to 
positive frequencies x only. Had we carried out those expansions over both positive and negative 
frequencies w , we would have obtained, instead of fl50|) and flSTl) . similar expressions 

2 u^uiu) = J{u)) o-^yiuj) and o-'yuioj) = 2 p,{ijj) u^uiuj) . (64) 

For positive cu, these would simply coincide with f l50l) and flHT]) . if we rename a; as x • For 
negative uj = —Xi the resulting expressions would read as 

2M7i.(-x) = J{.-x) o^-yvi.-X) and ^^v{,-x) = H - x) Uj^i - x) , (65) 

^ Churkin (1998) used functions which he caUed fc;(t) and which were, due to a difference in notations, the 
same as our ki (r) . 



14 



where we stick to the agreement that x always stands for a positive quantity. In accordance 
with ( l30l) . complex conjugation of ( l65|) would then return us to ( 164|) . 

Physically, the negative-frequency components of the stress or strain are nonexistent. If 
brought into consideration, they are obliged to obey (l30l) and, thus, should play no role, except 
for a harmless renormalisation of the Fourier components in (l32l) . 

When we say that the physically measurable stress (y^y{t) is equal to Xl"^^ (^-yuix)^^ 
it is unimportant to us whether the ^'Contribution in <J^^(t) comes from the term a'^;y(x) e ^ 

only, or also from the term a^^i ^x)^ ■ Indeed, the real part of the latter is a clone of the 
real part of the former (and it is only the former term that is physical). However, things remain 
that simple only for the stress and the strain. 

As we emphasised in subsection 13.41 the situation with the potentials is drastically different. 

While the physically measurable potential U{t) is still equal to ^TZe lf{x) e ^ , it is now 
important to distinguish whether the ^-contribution in U(t) comes from the term tJ^uix) ^ ^ 

or from the term U{ — x)g , ot perhaps from both. Although the negative mode — x would 
bring the same input as the positive mode x ? these inputs will contribute differently into the 



tidal torque. As can be seen from f l285p . the secular part of the tidal torque is proportional to 
sin , where = w^^p^ ^^impq ' with the time lag At^^^^ being positively defined - see formula 
(11091) . Thus the secular part of the tidal torque explicitly contains the sign of the tidal mode 

For this reason, as explained in subsection 13.41 a more accurate form of formula ( l63l) should 

be: 

f/,M = k{uj)W,{u) , (66) 

where u can be of any sign. 

If however, we pretend that the potentials depend on the physical frequency x = l^^l only, 
i.e., if we always write U{uj) as U{x) > then (163|) must be written as: 

Ui{.x) = kix) Wlix) , when X = I'^l for a; > , (67a) 

and 

Uiix) = k;{x) Wi{x) , when x = 1^1 for < . (67b) 

Unless we keep this detail in mind, we shall get a wrong sign for the Impq component of the 
torque after the despinning secondary crosses the appropriate commensurability. (We shall, of 
course, be able to mend this by simply inserting the sign sgau^^^^ by hand.) 



4.3 The complex Love number as a function of the complex compli- 
ance 

While the static Love numbers depend on the static rigidity modulus /i via ([3]), it is not readily 
apparent that the same relation interconnects ki{x) with /i(x) , the quantities that are the 
Fourier components of the time-derivatives of k2{t') and fi{t') . Fortunately, the correspondence 
principle (discussed in Appendix [D| tells us that, in many situations, the viscoelastic operational 
moduli fiix) or J(x) obey the same algebraic relations as the elastic parameters n or J . This 
is why, in these situations, the Fourier or Laplace transform of our viscoelastic equations will 
mimic (I228al - [228bp . except that all the functions will acquire overbars: a = 2 ft u , etc. 
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So their solution, too, will be Ui = kiWi , with ki retaining the same functional dependence on 
p, R, and fi as in ([3]), except that now /i will have an overbar: 



hix) 



1 3 1 



2(/-l) (2/^ + 4/ + 3)/i(x) 2(/ - 1) 1 + 

3 1 3 J(x) 



2(/ - 1) 1 + J/J{x) 2{l - 1) J(x) + AiJ 

Here the coefficients A; are defined via the unrelaxed quantities p = fi{0) = 1/J = 1/J{0) 
in the same manner as the static Ai were introduced through the static (relaxed) n = 1/J in 
formulae ([3]). 

The moral of the story is that, at low frequencies, each ki depends upon /2 (or upon J) in 
the same way as its static counterpart ki depends upon the static fj, (or upon the static J). 
This happens, because at low frequencies we neglect the acceleration term in the equation of 
motion fl231bp . so this equation still looks like ( ]228bp . 

Representing a complex Love number as 

kix) = ne[kix)] + ilrn[kix)] = |^Kx)l e""'^""^ (69) 
we can write for the phase lag ei{x) '■ 

= - 

ne [hix)] 



_ Xm [ki{x)\ 
tanei(x) = - ^ \j , A (70) 



or, equivalently: 

\ki{,x)\ sinei(x) = -Xm\ki{x)\ ■ (71) 

The products |fc/(x)| sinei(x) standing on the left-hand side in (ITTi) emerge also in the Fourier 
series for the tidal potential. Therefore it is these products (and not ki/Q ) that should enter 
the expansions for forces, torques, and the damping rate. This is the link between the body's 
rheology and the history of its spin: from J(x) to kiix) to |A;;(x)| sine(x) , the latter being 
employed in the theory of bodily tides. 

Through simple algebra, expressions ( 168|) entail: 

\ki[x)\ smei(x) = -Xm[ki[x)\ - 



2 (/ - 1) ( 7^e [J{x)] + Aijy + [Xm [J(x)] )' 

As we know from subsections 13.41 and 14. 2[ formulae ([7n]-[72]) should be used with care. Since 
in reality the potential U and therefore also ki are functions not of x but of u , then formulae 
(1721) should be equipped with multipliers sgn oj^^^^ , when plugged into the expression for the 
Impq component of the tidal force or torque. This prescription is equivalent to ([6] 



4.4 Should we write ki^pq and eimpq , or would ki and e/ be enough? 

In the preceding subsection, the static relation ([2]) was generalised to evolving settings as 

Uimpg{r,t) = {-\ ki{t)Wim„{R, r\t') , (73) 
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where Impq is a quadruple of integers employed to number a Fourier mode in the Darwin- 
Kaula expansion f llOOp of the tide, while Uimpq{r, t) and Wimpq{R, r*, t') are the harmonics 
containing cos{ximpqt - (^impq) and cos{ximpqt') correspondingly. 
One might be tempted to generalise ([2]) even further to 

Ulmpq{r, t) = l—j klmpq{t) Wlmpq{R , r*,t') , 

with the Love operator (and, consequently, its kernel, the Love function) bearing dependence 
upon m, p, and q . Accordingly, fl63|l would become 

Ulmpqix) = klmpqix) Wlmpq{x) ■ (74) 

Fortunately, insofar as the Correspondence Principle is valid, the functional form of the function 
klmpqix) depends upon / only and, thus, can be written down simply as kiiximpq) ■ know 
this from the considerations offered after equations fl228al - [228b|) . There we explained that ki 
depends on x = Ximpq only via J(x) , while the functional form of ki bears no dependence on 
X = Ximpq and, therefore, no dependence on m, p, q . 

The phase lag is often denoted as eimpq , a time-honoured tradition established by Kaula 
(1964). However, as the lag is expressed through ki via (!70l) . we see that all said above about 
ki apphes to the lag too: while the functional form of the dependency eimpqix) Kiay be different 
for different Is, it is invariant under the other three integers, so the notation eiiximpq) would be 
more adequate. 

It should be mentioned, though, that for bodies of pronounced non-sphericity coupling be- 
tween the spherical harmonics furnishes the Love numbers and lags whose expressions through 
the frequency, for a fixed I , have different functional forms for different m, p, q . In these cases, 
the notations kimpg and eimpq become necessary (Smith 1974; Wahr 1981a,b,c; Dehant 1987a,b). 
For a slightly non-spherical body, the Love numbers differ from the Love numbers of the spherical 
reference body by a term of the order of the flattening, so a small non-sphericity can usually be 
neglected. 



4.5 Rigidity vs self-gravitation 

For small bodies and small terrestrial planets, the values of vary from about unity to dozens 
to hundreds. For example, A2 is about 2 for the Earth (Efroimsky 2012), about 20 for Mars 
(Efroimsky & Lainey 2007), about 80 for the Moon (Efroimsky 2012), and about 200 for lapetus 
( Castillo- Rogez et al. 2011). For superearths, the values will be much smaller than unity, though. 
Insofar as 

Al > 1 , (75) 



one can approximate (l68ll with 



Ax) 



except in the closest vicinity of an Impq resonance, where the tidal frequency Ximpq approaches 
nil, and J diverges for some rheologies - like, for example, for those of Maxwell or Andrade. 
Whenever the approximate formula ( |76ll is applicable, we can rewrite ( |70l) as 



tane(x) - - ^H^^^ ^ - ^H^]^ = tan5(x) , (77) 

^-"^ ne[kiix)] ne[Jix)' ^ ' ^ ' 
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wherefrom we readily deduce that the phase lag e{x) of the tidal frequency x coincides with 
the phase lag of the complex compliance: 

e{x) ^ Six) , (78) 

provided x is ^^t too close to nil (i.e., provided we are not too close to the commensurability) . 
This way, insofar as the condition (71) is fulfilled, the component Ui{x) of the primary's potential 
lags behind the component Wi{x) of the perturbed potential by the same phase angle as the 
strain lags behind the stress at frequency x ^ sample of the material. Dependent upon the 
rheology, a vanishing tidal frequency may or may not limit the applicability of (71) and thus 
cause a considerable difference between e and 6 . 

In other words, the suggested approximation is valid insofar as changes of shape are deter- 
mined solely by the local material properties, and not by self-gravitation of the object as a whole. 
Whether this is so or not - depends upon the rheological model. For a Voigt or solid in 

the limit of x , we have J(x) — )■ J , so the zero- frequency limit of ki{x) is the static Love 
number ki = \k{0) \ . In this case, approximation (176]- 1781) remains applicable all the way down 
to X = . For the Maxwell and Andrade models, however, one obtains, for vanishing frequency: 
J{x) ~ ^/ivx) 5 whence ^ rjx and k2{x) approaches the hydrodynamical Love number 
kf"^"^ = ?>/2. 

We see that, for the Voigt and SAS models, approximation f lTHl) can work, for ^ 1, at 
all frequencies, because the condition ^ 1 can be set for all frequencies. For the Maxwell 
and Andrade solids, this condition holds only at frequencies larger than t~^A~^ = ^ ^^^^ 
so does the approximation (!78l) . Indeed, at frequencies below this threshold, self-gravitation 
"beats" the local material properties of the body, and the behaviour of the tidal lag deviates 
from that of the lag in a sample. This deviation will be indicated more clearly by formula (IMjl 
in the next section. The fact that, for some models, the tidal lag e deviates from the material 
lag angle 6 at the lowest frequencies should be kept in mind when one wants to explore crossing 
of a resonance. 

A standard caveat is in order, concerning formulae (1761 - 1781) . Since in reality the potential U 
is a function of u and not x ? our illegitimate use of x should be compensated by multiplying 
the function ^liXimpq) with sgncu^^p^ , when the lag shows up in the expression for the tidal 
force or torque. 

4.6 The case of inhomogeneous bodies 

Tidal dissipation within a multilayer near-spherical body is studied through expanding the in- 
volved fields over the spherical harmonics in each layer, setting the boundary conditions on the 
outer surface, and using the matching conditions on boundaries between layers. This formalism 
was developed by Alterman et al (1959). An updated discussion of the method can be found in 
Sabadini & Vermeersen (2004). For a brief review, see Legros et al (2006). 

Calculation of tidal dissipation in a Jovian planet is an even more formidable task (see Remus 
et al. 2012a and references therein). However dissipation in a giant planet with a solid core may 
turn out to be approachable by analytic means (Remus et al. 2011, 2012b). 

® The acronym SAS stands for the Standard Anelastic Solid, which is another name for the Hohenemser-Prager 
viscoelastic model. See the Appendix for details. 
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5 Dissipation at different frequencies 



5.1 The data collected on the Earth: in the lab, 

over seismological basins, and through geodetic measurements 

In Efroimsky & Lainey (2007), we considered the generic rheological model 

Q = i£xf , (r9a) 

where x is the tidal frequency and £^ is a parameter having the dimensions of time. The physical 
meaning of this parameter is elucidated in Ibid.. Under the special choice of a = —1 and for 
sufficiently large values of Q , this parameter coincides with the time lag At which, for this 
special rheology, turns out to be the same at all frequencies. 

Actual experiments register not the inverse quality factor but the phase lag between the 
reaction and the action. So the empirical law should rather be written down as 

{Sxf , (r9b) 



sin 5 



which is equivalent to fl79ap . provided the Q factor is defined there as Q energy ^^"^ '^'-'^ Qwork 
- see subsection 12.21 for details. 



The applicability realm of the empirical power law (179|) is remarkably broad - in terms of 
both the physical constituency of the bodies and their chemical composition. Most intriguing is 
the robust universality of the values taken by the index a for very different materials: between 
0.2 and 0.4 for ices and silicates, and between 0.14 and 0.2 for partial melts. Historically, two 
communities independently converged on this form of dependence. 



In the material sciences, the rheological model fl86|) . wherefrom the power law fl79b|) stems, 
traces its lineage to the groundbreaking work by Andrade (1910) who explored creep in metals. 
Through the subsequent century, this law was found to be applicable to a vast variety of other 
materials, including minerals (Weertman & Weertman 1975, Tan et al. 1997) and their partial 
melts (Fontaine et al. 2005). As recently discovered by McCarthy et al. (2007) and Castillo- 
Rogez (2009), the same law, with almost the same values of a , also applies to ices. The result 
is milestone, taken the physical and chemical differences between ices and silicates. It is agreed 
upon that in crystalline materials the Andrade regime can find its microscopic origin both in 
the dynamics of dislocations (Karato & Spetzler 1990) and in the grain-boundary diffusional 
creep (Gribb & Cooper 1998). As the same behaviour is inherent in metals, silicates, ices, and 
even glass-polyester composites (Nechada et al. 2005), it should stem from a single underlying 
phenomenon determined by some principles more general than specific material properties. An 
attempt to find such a universal mechanism was undertaken by Miguel et al. (2002). See also 
the theoretical considerations offered in Karato & Spetzler (1990). 

In seismology, the power law fl7^ became popular in the second part of the XX*'^ century, with 
the progress of precise measurements on large seismological basins (Mitchell 1995, Stachnik et al. 
2004, Shito et al. 2004). Further confirmation of this law came from geodetic experiments that 
included: (a) satellite laser ranging (SLR) measurements of tidal variations in the J2 component 
of the gravity field of the Earth; (b) space-based observations of tidal variations in the Earth's 
rotation rate; and (c) space-based measurements of the Chandler Wobble period and damping 
(Benjamin et al. 2006, Eanes & Bettadpur 1996, Eanes 1995). Not surprisingly, the Andrade 
law became a key element in the recent attempt to construct a universal rheological model of 
the Earth's mantle (Birger 2007). This law also became a component of the non-hydrostatic- 
equilibrium model for the zonal tides in an inelastic Earth by Defraigne & Smits (1999), a 
model that became the basis for the lERS Conventions (Petit & Luzum 2010). While the lab 
experiments give for a values within 0.2 — 0.4 , the geodetic techniques favour the interval 
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0.14 — 0.2. This minor discrepancy may have emerged due to the presence of partial melt in 
the mantle and, possibly, due to nonlinearity at high bounding pressures in the lower mantle. 
The universality of the Andrade law compels us to assume that (!79|) works equally well for other 
terrestrial bodies. Similarly, the applicability of ( !79|) to samples of ices in the lab is likely to 
indicate that this law can be employed for description of an icy moon as a whole. 

Karato & Spetzler (1990) argue that at frequencies below a certain threshold xo anelasticity 
gives way to purely viscoelastic behaviour, so the parameter a becomes close to unityj^ For 
the Earth's mantle, the threshold corresponds to the time-scale about a year or slightly longer. 
Although in Karato & Spetzler (1990) the rheological law is written in terms of 1/Q , we shall 
substitute it with a law more appropriate to the studies of tides: 

ki sinei = {S x) ? where p = 0.2 — 0.4 for x > Xo 

and p ~ 1 for X < Xo , (80) 

X being the frequency, and xo being the frequency threshold below which viscosity takes over 
anelasticity. 

The reason why we write the power scaling law as (IHOl) and not as (!79|) is that at the lowest 
frequencies the geodetic measurements give us actually ki sin e/ = — Zm [ k^ (x) ] and not the lag 
angle 5 in a sample (e.g., Benjamin et al. 2006). For this same reason, we denoted the exponents 
in ( |79ll and (IHOll with different letters, a and p . Below we shall see that these exponents do not 
always coincide. Another reason for giving preference to flHOl) is that not only the sine of the lag 
but also the absolute value of the Love number is frequency dependent. 

5.2 Tidal damping in the Moon, from laser ranging 

Fitting of the LLR data to the power scaling law (1791) . which was carried out by Williams et 
al. (2001), has demonstrated that the lunar mantle possesses quite an abnormal value of the 
exponent: —0.19. A later reexamination in Williams et al. (2008) rendered a less embar- 
rassing value, — 0.09 , which nevertheless was still negative and thus seemed to contradict our 
knowledge about microphysical damping mechanisms in minerals. Thereupon, Williams & Boggs 
(2009) commented: 

" There is a weak dependence of tidal specific dissipation Q on period. The Q increases from 
~ 30 at a month to ~ 35 at one year. Q for rock is expected to have a weak dependence on tidal 
period, but it is expected to decrease with period rather than increase. The frequency dependence 
of Q deserves further attention and should be improved." 

While there always remains a possibility of the raw data being insufficient or of the fitting 
procedure being imperfect, the fact is that the negative exponent obtained in Ibid, does not 
necessarily contradict the scaling law (1791) proven for minerals and partial melts. Indeed, the 
exponent obtained by the LLR Team was not the a from (1791) but was the p from ( IHOil . The 
distinction is critical due to the difference in frequency-dependence of the seismic and tidal dissi- 
pation. It turns out that the near- viscous value p ~ 1 from the second line of flHOl) . appropriate 
for low frequencies, does not retain its value all the way to the zero frequency. Specifically, in 
subsection 15.41 we shall see that at the frequency — ^t- (where r^^ = r^/yU is the Maxwell 

time, with rj and fi being the lunar mantle's viscosity and rigidity), the exponent p begins 

^ This circumstance was ignored by Defraigne & Smits (1999). Accordingly, if the claims by Karato & Spetzler 
(1990) are correct, the table of corrections for the tidal variations in the Earth's rotation in the lERS Conventions 
is likely to contain increasing errors for periods of about a year and longer. 

This detail is missing in the theory of the Chandler wobble of Mars, by Zharkov & Gudkova (2009). 
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to decrease with the decrease of the frequency. As the frequency becomes lower, p changes its 
sign and eventually becomes — 1 in a close vicinity of x = • This behaviour follows from 
calculations based on a realistic rheology (see formulae (p2l - below), and it goes along well 
with the evident physical fact that the average tidal torque must vanish in a resonances In 
subsection 15.71 comparison of this behaviour with the LLR results will yield us an estimate for 
the mean lunar viscosity. 

5.3 The Andrade model as an example of viscoelastic behaviour 

The complex compliance of a Maxwell material contains a term J = J(0) responsible for the 
elastic part of the deformation and a term — describing the viscosity. Whatever other terms 
get incorporated into the compliance, these will correspond to other forms of hereditary reaction. 
The available geophysical data strongly favour a particular extension of the Maxwell approach, 
the Andrade model (Cottrell & Aytekin 1947, Duval 1976). In modern notations, the model can 
be expressed as 

j(t-t') = [J + l3{t-tT + V'Ht~t')] ^{t-t') , (81) 

a being a dimensionless parameter, /3 being a dimensional parameter, t] denoting the steady- 
state viscosity, and J standing for the unrelaxed compliance, which is inverse to the unrelaxed 
rigidity: J = J(0) = l//i(0) = l//x. We see that (ISTl) is the Maxwell model amended with an 
extra term of a hereditary nature. 

A simple example illustrating how the model works is rendered by deformation under constant 
loading. In this case, the anelastic term dominates at short times, the strain thus being a convex 
function of t (the so-called primary or transient creep). As time goes on and the applied loading 
is kept constant, the viscous term becomes larger, and the strain becomes almost linear in time 
- a phenomenon called the secondary creep. 

Remarkably, for all minerals (including ices) the values of a belong to the interval from 0.14 
through 0.4 (more often, through 0.3) - see the references in subsection 15.11 above. The other 
parameter, /3 , may be rewritten as 

P = J r;" = f,-' r;° , (82) 

the quantity having dimensions of time. This quantity is the timescale associated with the 
Andrade creep, and it may be termed as the "Andrade time" or the "anelastic time" . It is clear 
from (1821) that a short makes the anelasticity more pronounced, while a long makes the 
anelasticity weakQ 

It is known from Castillo-Rogez et al. (2011) and Castillo-Rogez & Choukroun (2010) that 
for some minerals, within some frequency bands, the Andrade time gets very close to the Maxwell 
time: 

^ r,, =^ P ^ J r7 = J^"" r/"" = /.""^ r/"" , (83) 

^'^ For example, the principal tidal torque t^^^^^ = ''"2200 ^-cting on a secondary must vanish when the secondary 
is crossing the synchronous orbit. Naturally, this happens because p becomes — 1 in the close vicinity of 

X2200 ~ ^ ■ 

^^As long as we agree to integrate over t — t' e [ 0, 00 ) , the terms /3{t~t')" and rj^^ {t — t') can do without 
the Heaviside step-function Q{t — t'). We remind though that the first term, J, does need this multiplier, so 
that insertion of ([8T|) into (|43|) renders the desired J 5{t — t') under the integral, after the differentiation in (j43|) 
is performed. 

While the Andrade creep is likely to be caused by "unpinning" of jammed dislocations (Karato & Spetzler 
1990, Miguel et al 2002), it is not apparently clear if the Andrade time can be identified with the typical time of 
unpinning of defects. 
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where the relaxation Maxwell time is given by: 

_ V 



Tj J 



On general grounds, though, one cannot expect the anelastic timescale and the viscoelastic 
timescale r^^ to coincide in all situations. This is especially so due to the fact that both these 
times may possess some degree of frequency-dependence. Specifically, there exist indications 
that in the Earth's mantle the role of anelasticity (compared to viscoelasticity) undergoes a 
decrease when the frequencies become lower than 1/yr - see the miscrophysical model suggested 
in subsection 5.2.3 of Karato & Spetzler (1990). It should be remembered, though, that the 
relation between and r^^ may depend also upon the intensity of loading, i.e., upon the 
damping mechanisms involved. The microphysical model considered in Ihid. was applicable 
to strong deformations, with anelastic dissipation being dominated by dislocations unpinning. 
Accordingly, the dominance of viscosity over anelasticity ( ^ r^^ ) at low frequencies may be 
regarded proven for strong deformations only. At low stresses, when the grain-boundary diffusion 
mechanism is dominant, the values of and t^^^ may remain comparable at low frequencies. 
The topic needs further research. 

In terms of the Andrade and Maxwell times, the compliance becomes: 



J{t-t') = J 



1 + 



t -t' 



t-t' 



Q{t-t') 



15) 



In the frequency domain, compliance (18511 will look: 



Hx) = J + /3(^x)-"r(i + «) 

= J [1 + (^xrj-°r(l + a) - i{xr,,)-^] 



f86a) 



(86b) 



X being the frequency, and F denoting the Gamma function. The imaginary and real parts of 
the complex compliance are: 



lm[J{x)] 



— - sin(^) T{a + 1) 

VX \ 2 J 

■HxrJ-' - J(xrJ- sin(^) r(a + 1) 



and 



^e[J(x)] 



J + X cos 



a 71 



r(a + 1) 



- cos ( ^ ) Via + 1) 



= J + J{XT^ 

whence we obtain the following dependence of the phase lag upon the frequency: 

13 sin ( ^) r (a + 



tanS{x) 
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Here z is the dimensionless frequency defined as 

z = X^j, = X^^M ^ ) (90) 
while C is a dimensionless parameter of the Andrade model: 

C = ^ . (91) 

5.4 Tidal response of viscoelastic near-spherical bodies obeying the 
Andrade and Maxwell models 

An Im-pq term in the expansion for the tidal torque is proportional to the factor (x) sin (x) = 
^iOiim'pq)\ ^^"^^i^Ximpq) ' Heuce the tidal response of a body is determined by the frequency- 
dependence of these factors. 

Combining fl72p with (186 p , and keeping in mind that A; 3> 1 , it is easy to write down the 
frequency- dependencies of the products |A;j(x)| sinej(x) • Referring the reader to Appendix IE.2I 
for details, we present the results, without the sign multiplier. 

• In the high-frequency band: 
Ux)\ sine.ix) - sin(^) r(a + 1) {r,.^ xT" , for X » r"^ 

For small bodies and small terrestrial planets (i.e., for ^ 1 ), the boundary between the high 
and intermediate frequencies turns out to be 



_ -1 

Xhi '^m ^ 

For large terrestrial planets (i.e., for ^ 1 ) the boundary frequency is 

A.HI 'a M ^ 

At high frequencies, anelasticity dominates. So, dependent upon the microphysics of the mantle, 
the parameter C, may be of order unity or slightly lower. We say slightly, because we expect 
both anelasticity and viscosity to be present near the transitional zone. (A too low C, would 
eliminate viscosity from the picture completely.) This said, we may assume that the boundary 
Xhi is comparable to r^^^ for both small and large solid objects. This is why in (|92|) we set the 
inequality simply as X ^ Tj^^ • 

• In the intermediate-frequency band: 
|fc,(x)| sine,(x) ^ JlfZY) [A^llf ' ^m' » ^ » ^m' (A + 1)"' • 



While the consideration in the Appendix IE. 21 renders C ^ " for the upper bound, here we 
approximate it with in understanding that C, does not differ from unity too much near the 
transitional zone. Further advances of rheology may challenge this convenient simplification. 
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for (A + 1) ' > X 



Scaling laws ( l92l) and ( |93l) mimic, up to constant factors, the frequency-dependencies of 
|J(x)| sin5(x) = — Xm[J(x)] at high and low frequencies, correspondingly, - this can be seen 
from (1871). 

Expression flM|) however shows a remarkable phenomenon inherent only in the tidal lagging, 
and not in the lagging in a sample of material: at frequencies below t~^{A^ + 1)~'^ = ^(A^ + l)"^, 

the product smei{x) changes its behaviour and becomes linear in x- 

While elsewhere the sine;(x) factor increases with decreasing x, it changes its 

behaviour drastically on close approach to the zero frequency. Having reached a finite maximum 
at about x = '^Z^i^i + ^^e said factor begins to scale linearly in x as x approaches zero. 
This way, the factor sinej(x) decreases continuously on close approach to a resonance, 

becomes nil together with the frequency at the point of resonance. So neither the tidal torque 
nor the tidal force explodes in resonances. In a somewhat heuristic manner, this change in the 
frequency- dependence was pointed out, for / = 2 , in Section 9 of Efroimsky & Williams (2009). 

5.5 Example 

Figured] shows the absolute value, ^2 = |^2(x)l > well as the real part, TZe [k2{x)] = ^2 cos €2 , 
and the negative imaginary part, — Xm [k2ix)] = ^2 sin €2 , of the complex quadrupole Love 
number. Each of the three quantities is represented by its decadic logarithm as a function of 
the decadic logarithm of the forcing frequency x (given in Hz). The curves were obtained by 
insertion of formulae (IHTl - IHHl) into ( |68l) . As an example, the case of —Xm [^2(x)] is worked 
out in Appendix IE. 21 see formulae (12521 - [HHI) • 

Both in the high- and low-frequency limits, the negative imaginary part of k2{x) 1 given on 
Figure [1] by the red curve, approaches zero. Accordingly, over the low- and high-frequency bands 
the real part (the green line) virtually coincides with the absolute value (the blue line). 

While on the left and on the close right of the peak, dissipation is mainly due to viscosity, 
friction at higher frequencies is mainly due to anelasticity. This switch corresponds to the change 
of the slope of the red curve at high frequencies (for our choice of parameters, at around 10~^ 
Hz). This change of the slope is often called the elbow. 

Figured] was generated for A2 = 80.5 and t^^j = 3.75 x 10^ s. The value of A2 corresponds 
to the Moon modeled by a homogeneous sphere of rigidity /i = 0.8 x 10^^ Pa. Our choice of the 
value of =1]/ n corresponds to a homogeneous Moon with the said value of rigidity and with 
viscosity set to be = 3 x 10^^ Pa s. The reason why we consider an example with such a low 
value of 7] will be explained in subsection 15. 7[ Finally, it was assumed for simplicity that C = 1 , 
i.e., that = Tj^ . Although unphysical at low frequencies, this simplification only slightly 
changes the shape of the "elbow" and exerts virtually no influence upon the maximum of the red 
curve, provided the maximum is located well into the viscosity zone. 

5.6 Crossing a resonance — with a chance for entrapment 

As ever, we recall that in the expansion for the tidal torque factors ^3% - [Ml) should appear in 
the company of multipliers sgn uj . For example, the factor describing dissipation near an 
Impq resonance will enter the expansions as 

l^i(X«mpg)l ^i^^«(Ximpg) ^?>^^lmpq ~ 2 (/ _ 1) ^Irnpq ^S^^lmpq ~ 2 {I — 1) ^'^^ ^Impq " 
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• In the low-frequency band: 
3 

\kiix)\ sineiix) ^ 2(1-1) ^ 




Figure 1: Tidal response of a homogeneous spherical Andrade body, set against the decadic 
logarithm of the forcing frequency x iy^ Hz). The blue curve renders the decadic logarithm 
of the absolute value of the quadrupole complex Love number, lg/c2 = lg|^2(x)l • The green 
and red curves depict the logarithms of the real and the negative imaginary parts of the Love 
number: IgT^e [A;2(x)] = lg(fc2Cose2) and Ig {— Xm [A;2(x)] } = Ig (— ^2 sin €2) , accordingly. 
The change in the slope of the red curve (the "elbow"), which takes place to the right of the 
maximum, corresponds to the switch from viscosity dominance at lower frequencies to anelasticity 
dominance at higher frequencies. The parameters A2 and t^^^ were given values appropriate 
to a homogeneous Moon with a low viscosity, as described in subsection 15.71 The plots were 
generated for an Andrade body with = 1 at all frequencies. Setting the body Maxwell at lower 
frequencies will only slightly change the shape of the "elbow" and will have virtually no effect 
on the maximum. 
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Naturally, the Impq term of the torque depends on ujf^p^ and not just on Ximpg ~ l^impgl ■ '^^^^ 
term then goes continuously through zero, and changes its sign as the Impq resonance is crossed 
(i.e., as oJi^pg goes through nil and changes its sign). 

Formula (!95|) tells us that an Impq component of the tidal torque goes continuously through 
zero as the satellite is traversing the commensurability which corresponds to vanishing of a tidal 
frequency Ximpq ■ This gets along with the physically evident fact that the principal (i.e., 2200 ) 
term of the tidal torque should vanish as the secondary approaches the synchronous orbit. 

It is important that a Impq term of the torque changes its sign and thus creates a chance for 
entrapment. As the value of an Impq term of the torque is much lower than that of the principal, 
2200 term, we see that a perfectly spherical body will never get stuck in a resonance other than 
2200. (The latter is, of course, the 1 : 1 resonance, i.e., the one in which the principal term 
of the torque vanishes.) However, the presence of the triaxiality-generated torque is known to 
contribute to the probabilities of entrapment into other resonances (provided the eccentricity is 
not zero). Typically, in the literature they consider a superposition of the triaxiality-generated 
torque with the principal tidal term. We would point out that the "trap" shape of the Impq 
term (!95|) makes this term relevant for the study of entrapment in the Impq resonance. In some 
situations, one has to take into account also the non-principal terms of the tidal torque. 



5.7 Comparison with the LLR results 

As we mentioned above, fitting of the lunar laser ranging (LLR) data to the power law has resulted 
in a very small negative exponent p = —0.19 (Williams et al. 2001). Since the measurements 
of the lunar damping described in Ibid, rendered information on the tidal and not seismic 
dissipation, those results can and should be compared to the scaling law ([92] -[91]). As the small 
negative exponent was devised from observations over periods of a month to a year, it is natural 
to presume that the appropriate frequencies were close to or slightly below the frequency — ^t— 

at which the factor k2 sin 62 has its peak: 



3xl0'.«.0.1yr = r„^,= ^^,= , (96) 

as on Figure [H Hence, were the Moon a uniform viscoelastic body, its viscosity would be only 

ri = 3 X 10^^ Pa s . (97) 

For the actual Moon, the estimate means that the lower mantle contains a high percentage of 
partial melt, a fact which goes along well with the model suggested in Weber et al. (2011), and 
which was anticipated yet in Williams et al. (2001) and Williams & Boggs (2009), following an 
earlier work by Nakamura et al. (1974). 



6 The polar component of the tidal torque acting on the 
primary 

Let vector r = (r. A, 0) point from the centre of the primary toward a point-like secondary 
of mass Msec ■ Associating the coordinate system with the primary, we reckon the latitude 
from the equator. Setting the coordinate system to corotate, we determine the longitude A 
from a fixed meridian. The tidally induced component of the primary's potential, U , can be 
generated either by this secondary itself or by some other secondary of mass M*^^ located at 
r* = (r*. A*, 0*) . In either situation, the tidally induced potential U generates a tidal force 
and a tidal torque wherewith the secondary of mass M^ec acts on the primary. 
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The scope of this paper is hmited to low values of i. When the role of the primary is played 
by a planet, the secondary being its satellite, i is the satellite's inclination. When the role of the 
primary is played by the satellite, the planet acting as its secondary, i acquires the meaning of 
the satellite's obliquity. Similarly, when the planet is regarded as a primary and its host star is 
treated as its secondary, i is the obliquity of the planet. In all these cases, the smallness of i 
indicates that the tidal torque acting on the primary can be identified with its polar component, 
the one orthogonal to the equator of the primary. The other components of the torque will be 
neglected in this approximation. 

The polar component of the torque acting on the primary is the negative of the partial 
derivative of the tidal potential, with respect to the primary's sidereal angle: 

r(r^ = - M,,, , (98) 

6 standing for the primary's sidereal angle. This formula is convenient when the tidal potential 
U is expressed through the secondary's orbital elements and the primary's sidereal ang 

Here and hereafter we are deliberately referring to a primary and a secondary in lieu of 
a planet and a satellite. The preference stems from our intention to extend the formalism to 
setting where a moon is playing the role of a tidally-perturbed primary, the planet being its 
tide-producing secondary. Similarly, when addressing the rotation of Mercury, we interpret the 
Sun as a secondary that is causing a tide on the effectively primary body. Mercury. 



7 The tidal potential 

7.1 Darwin (1879) and Kaula (1964) 

The potential produced at point i? = (i? , A , 0) by a secondary body of mass M*, located at r* = 
(r*. A* , (p*) with r* > R, is given by ([T]). When a tide-raising secondary located at r* distorts 
the shape of the primary, the potential generated by the primary at some exterior point r gets 
changed. In the linear approximation, its variation is given by ([2]). Insertion of ([1]) into ([2]) entails 

U{r) = -GMl^^Y^h ^(2-MP,^(sin0)P,^(sinr) cosm(A-A*) . 

1=2 r r* ^=oy^ + ^)- 

A different expression for the tidal potential was offered by Kaula (1961, 1964), who developed 
a powerful technique that enabled him to switch from the spherical coordinates to the Kepler 
elements (a*, e*, i*, fl*, u*, A^* ) and (a, e, i, fl, u, J^) of the secondaries located at r* and 

13 Were the potential written down in the spherical coordinates associated with the primary's equator and 
corotating with the primary, the polar component of the tidal torque could be calculated with aid of the expression 

r(f) = M... ^ 

derived, for example, in Williams & Efroimsky (2012). That the expression agrees with can be seen from the 
formula 

X = -e + n + Lu + i^ + 0{i^) = -e + n + Lj + M+2e smM + 0{e^) + 0{i^) , 

e , i ^ uj , , 1^ and A4 being the eccentricity, inclination, argument of the pericentre, longitude of the node, true 
anomaly, and mean anomaly of the tide-raising secondary. 
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r . Application of this technique to ( l99l) resuhs in 



(100) 



. a / a* \ a* J ^ (/ + m)! 
=2 ^ ^ ^ ^ m=0 ^ ^ 



Z oo Z oo 

p=0 g=— oo /i=0 j=— oo 



where 



v\^^^ = {I- 2p)uj* + (/ - 2p + q)M* + mn* , (101) 

= {l-2h)u + {l-2h + j)M + niQ , (102) 

9 = 9* being the sidereal angle, Gipg{e) signifying the eccentricity functions^ and Fimp{i) 
denoting the inclination functions (Gooding & Wagner 2008). 

Being equivalent for a planet with an instant response of the shape, fl99l) and fllOOp disagree 
when dissipation-caused delays come into play. Kaula's expression fllOOp . as well as its truncated, 
Darwin's version^ is capable of accommodating separate phase lags for each mode: 

1=2 ^ ^ ^ m=0 ^ ' 



(103) 



Spm ) Fimpji*) Gipq{e*)'^Fimh{i) ^ Gihj{e) COS [[v*^p^ - m9*) - { 

Vlmhj Tn9) (^Impql 

p=0 q=—oo h=0 j=—oo 

where 

eimpq = {I - 2p) U* + {I -2p + q) M* + m{Q* - 9*) Atimpq = OJ^rnpq ^Umpq = ± X*mpq ^Umpq (104) 

is the phase lag. The tidal mode w*^^^ introduced in fll04p is 

uj*^^^ ^ {l-2p)6j* + {l~2p + q)M* + m{tl* - 9*) , (105) 

while the positively-defined quantity 

X]mpq = I^LpJ = l(^-2p)u;* + {l-2p + q)M* + mih* - 9*) \ (106) 

is the actual physical Impq frequency excited by the tide in the primary. The correspond- 
ing positively-defined time delay Atimpq = Ati{ximpq) depends on this physical frequency, the 
functional forms of this dependence being different for different materials. 



Functions Gihj{e) coincide with the Hansen polynomials -^|;_2p^^) ^\^) ■ In Appendix [G| we provide a 

table of the Gihj (e) required for expansion of tides up to , inclusively. 

While the treatment by Kaula (1964) entails the infinite Fourier series (|100p . the development by Darwin 
(1879, 1880) renders its partial sum with |^|, \q\, \j\ < 2. For a simple introduction into Darwin's method see 
Ferraz-Mello et al. (2008). Be mindful that in Ibid, the convention on the notations f and r* is opposite to ours. 
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In neglect of the apsidal and nodal precessions, and also of A4o , the above formulae become: 

uimpg = {I -2p + q)n - m9 , (107) 

Xlmpq = I Wjmpg \ = \ {I - 2p + q) u - m 9 \ , (108) 

and 

eimpq = (^Impq^timpq = {I - 2p + q) 71 - TU 6 Atlmpq (109a) 

= Xlmpq AtliXlmpq) Sgn {I - 2p + q) U - 111 9 , (109b) 

Formulae f llOOp and fll03p constitute the pivotal result of Kaula's theory of tides. Importantly, 
Kaula's theory imposes no a priori constraint on the form of frequency-dependence of the lags. 

8 The Darwin torque 

As explained in Williams & Efroimsky (2012), the empirical model by MacDonald (1964), called 
MacDonald torque^ tacitly sets an unphysical rheology of the satellite's material. The rheology is 
given by (1791) with a = — 1 . More realistic is the dissipation law fl80p . An even more accurate 
and practical formulation of the damping law, stemming from the Andrade formula for the 
compliance, is rendered by (1921 - [9^ . These formulae should be inserted into the Darwin-Kaula 
theory of tides. 

8.1 The secular and oscillating parts of the Darwin torque 
8.1.1 The general formula 

Direct differentiation of (I103P with respect to — 6' will result in the expression^ 

1=2 ^ ^ ^ m=0 ^ ' p=0 q=-oo 

I oo 

y^^Fimhji) Gihj{e)ki sin 1^ W^mpg Vlfnfij (^Impq 

(110) 

h=0 j=—oo 

If the tidally-perturbed and tide-raising secondaries are the same body, then Msec = M*^^ , and 
all the elements coincide with their counterparts with an asterisk. Hence the differences 

* 

'^Impq " '^Imhj — 

{I - 2p + q) M* - {I - 2h + j)M + m{Q* -Q) + I {co* - u) - 2puj* + 2hu (111) 
get simphfied to 

vLpq - vimhj = {2h-2p + q-j)M* + {2h- 2p)uj* , (112) 

an expression containing both short-period contributions proportional to the mean anomaly, and 
long-period contributions proportional to the argument of the pericentre. 
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For justification of tliis operation, see Section 6 in Efroimsky & Williams (2009). 
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8.1.2 The secular, the purely-short-period, and the mixed-period parts of the 
torque 

Now we see that the terms entering series (IllOp can be spUt into three groups: 

(1) The terms, in which indices {p , q) coincide with [h , j) , constitute a secular part of 
the tidal torque, because in such terms vimhj cancel with f . This Ai- and w-independent 
part is furnished by 

OO -r-,2i + 1 i , , > , I OO 

T = ^2 G M,^^ , 777 ^Fi^p^i) Gip^{e) h sine^^pg . (113) 

1=2 O' m=0 ^ '' p=0 q=-oo 

(2) The terms with p = h and q ^ j constitute a purely short-period part of the torque: 

OO 2(+l / / , \ , I OO OO 
^ 2 j-^ I ^ 777 ) I ^ 2 ^ ^ 

7" = - 5^ 2 G M,,, Y ^ ^impi^) Yl Yl ^ipi(^) '^'pi(e) ki sin [ (g 

1=2 O rn=0 ^ '' p=0 q=-oo j=-oo 

-j) M - eimpq] . (114) 

(3) The remaining terms, ones with p ^ h , make a mixed-period part comprised of both short- 
and long-period terms: 

mixed 

T 

OO ^ J^^^^ ' f/ ' ' OO OO 

- ^2GMl^^^L^^m^F,„p(0$^F,^/.(z) E Gi,,{e) Gi,j{e) k, sin [ (2 /7 

1=2 m=0 ^ p=0 h=0 q=-OQ j=-oo 

h^p 

-2p + q - j) M* + i2h - 2p)L0* - eimpq] . (115) 

8.1.3 The / = 2 and / = 3 terms in the 0(7^) approximation 

For 1 = 2, index m will take the values 0,1,2 only. Although the rTi = terms enter the 
potential, they add nothing to the torque, because differentiation of f ll03p with respect to —9 
furnishes the m multiplier in flllOp . To examine the remaining terms, we should consider the 
inclination functions with subscripts (/777p) = (220) , (210) , (211) only: 

F22o(z) = 3 + 0(z2) , F2io(^) = ^ sinz + 0(7^) , ^211(7) = - ^ sin 7 + 0(7^) , (116) 

all the other F2mp{i) being of order O(z^) or higher. Thence for p = h (i.e., both in the secular 
and purely-short-period parts) it is sufficient, in the O(z^) approximation, to keep only the terms 
with ^2^20(0 5 ig^o^i^g those with -Ffio(0 -^2^11 (0 ■ We see that in the O(z^) approximation 

• among the I = 2 terms, both in the secular and purely short-period parts, 
only the terms with {Imp) = (220) are relevant. 

In the case of p ^ h , i.e., in the mixed-period part, the terms of the leading order in inclination 
are: Fimp{i) Fimh{i) = -p2io(«) -^2ii(0 and Fimp{i) Fimh{i) = -F2ii(z) F2io(z) , which happen 
to be equal to one another, and to be of order O(z^) . This way, in the O(i^) approximation. 
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the mixed-period part of the / = 2 component may be omitted. 



The inchnation functions Fi 



Imp 



310 ) -f^312 ; -^313 , -^320 , -^321 , -^322 ; -^323 ; -^^331 ; -^^332 ; -^333 



are of order 0{i) or higher. The terms containing the squares or cross-products of the these 
functions may thus be dropped. Specifically, the smallness of the cross-terms tells us that 

• the mixed-period part of the / = 3 component may be omitted. 

What remains is the terms containing the squares of functions 



In other words, 



2 + o(e) 



and 



^330 (0 



15 + Oii' 



;il7) 



• among the I = 3 terms, both in the secular and purely short-period parts, 
only the terms with {Imp) = (311) and {Imp) = (330) are important. 

All in all, for / = 2 and I = 3 the mixed-period parts of the torque may be neglected, in the 
O(i^) approximation. The surviving terms of the secular and the purely short-period parts will 
be developed up to , inclusively. 

8.2 Approximation for the secular and short-period parts of the tidal 
torque 

As we just saw, both the secular and short-period parts of the torque may be approximated with 
the following degree-2 and degree-3 components: 



r = r + r + o{e{R/af) 

1=2 1=3 V V / W 



r + 

(«mp)=(220) 



r + r 

(«mp)=(311) («mp)=(330)_ 



+ 0{ei^) + 0(e(i?/a)9) 



:il8) 



and 



f = + + 0{e{R/af) 



T + 

(Zmp)=(220) 



r + r 

(imp)=(311) (Zmp)=(330) _ 



0{ei') + 0{e{R/af) 



(119) 



were the / = 2 and / = 3 inputs are of the order {R/a)^ and {R/a)"^ , accordingly; while the 
/ = 4, 5, . . . inputs constitute O { e {R/a)^ ) . 

Expressions for T, , , , , T,, , , , , and T,, , , , are furnished bv formulae (12831) . 

^ («mp)=(220) ' (/mp)=(311) ' («mp)=(310) ^ ^' 
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(I285p . and ( I287p in Appendix [Hi As an example, here we provide one of these components: 
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+ 0(e«e) + O(z^e) 



(120) 



Here each term changes its sign on crossing the appropriate resonance. The change of the sign 
takes place smoothly, as the value of the term goes through zero - this can be seen from formula 
(195|) and from the fact that the tidal mode cu^^^^ vanishes in the Impq resonance. 

Expressions for 7j^^^)^(220) ' Vrt=(3ii) ' Vp)=(330) are given by formulae (1281- EUD 
in Appendix [H Although the average of the short-period part of the torque vanishes, this part 
does contribute to dissipation. Oscillating torques contribute also to variations of the surface of 
the tidally-distorted primary, the latter fact being of importance in laser-ranging experiments. 

Whether the short-period torque may or may not influence also the process of entrapment is 
worth exploring numerically. The reason why this issue is raised is that the frequencies n{q — j) 
of the components of the short-period torque are integers of n and thus are commensurate with 
the spin rate 9 near an A/B resonance, A and B being integer. It may be especially interesting 
to check the role of this torque when q — j = 1 and A/B = N is integer. 

The hypothetical role of the short-period torque in the entrapment and libration dynamics 
has never been discussed so far, as the previous studies employed expressions for the tidal torque, 
which were obtained through averaging over the period of the secondary's orbital motion. 



8.3 Librations 

Consider a tidally-perturbed body caught in a A : B resonance with its tide-raising companion, 
A and B being integer. Then the spin rate of the body is 



A 

9 = —n + %Ij , 

where the physical-libration angle is 



'1211 



(122) 
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oOp^ being the physical-libration frequency. The oscillating tidal torque exerted on the body is 
comprised of the modes 

^impq = {i-'^P + q)n-m6= ^-2p + q-^ n - m ip . (123) 

In those Impq terms for which the combination / — 2p + g — ^ m is not zero, the small 

quantity —mi/j may be neg lectedllll The remaining terms will pertain to the geometric libration. 
The phase lags will be given by the standard formula e^^^^ = uii^pg ^^i^pq ■ 

In those Impq terms for which the combination I — 2p + q — "^m vanishes, the physical- 



libration input ~mip is the only one left. Accordingly, the multiplier sin yv*^^^ — m9* j — 

^ "^impg — raO^ in the Impq term of the tidal torque will reduce to sin [ — m {ip* — if))] ^ 

—mijjAt = mipoOJpL At cosw^^t . Here the time lag At is the one corresponding to the 
physical-libration frequency which may be very different from the usual tidal frequencies 
for nonsynchronous rotation - see Williams & Efroimsky (2012) for a comprehensive discussion. 



9 Marking the minefield 

The afore-presented expressions for the secular and purely short-period parts of the tidal torque 
look cumbersome when compared to the compact and elegant formulae employed in the literature 
hitherto. It will therefore be important to explain why those simplifications are impractical. 



9.1 Perils of the conventional simplification 

Insofar as the quality factor is large and the lag is small (i.e., insofar as sine can be approximated 
with e ) , our expression (I282ap assumes a simpler form: 



q=~3 

where the error 0{e^) originates from sine = e + O(e^) . 

The simplification conventionally used in the literature ignores the frequency-dependence of 
the Love number and attributes the overall frequency-dependence to the lag. It also ignores the 
difference between the tidal lag e and the lag in the material, 6 . This way, the conventional 
simplification makes e obey the scaling law (179b|) . At this point, most authors also set a = —1 . 
Here we shall explore this approach, though shall keep a arbitrary. From the formula 

At, = £ X (125) 

Impq Impq ^ ' 

derived by Efroimsky & Lainey (2007) in the said approach, we see that the time lags are related 
to the principal-frequency lag ^^2200 



a+l 

At, = At ( ) . (126) 

Impq 2200 1 - 1 ^ ' 



Y 

Impq 



The physical-libration input —ma may be neglected in the expression for ijJi^p^ even when the magnitude 
of the physical libration is comparable to that of the geometric libration (as in the case of Phobos). 

^® Let .j^^ = (£%)" , where 8 is an empirical parameter of the dimensions of time, while e is small enough, 
so sine w e . In combination with e, = uj. At, and y, = \lo, | , these formulae entail (11251) . 

impq Irapq Irapq ^Impq ' tmpq ' ^ ^ 
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When the despinning is still going on and 6 ^ n , the corresponding phase lags are: 

Q + l 

e = At uj = - At Y I 1 = - e ^ "^^^"^ 

Impq Irnpq Impq 2200 ^Impq \ -v / 2200 



X, 



X, 



(127) 



Impq / \ Impq 

which helps us to cast the secular part of the torque into the following convenient form: 
(Q>10) . 
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O(zVg) + 0(g^3) + 0[ae^Q-^nlQ) (128a) 
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where the overall factor reads as: 
Z 
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+ M„ 



3nM,,^ k2 
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^2200 



Mprim and Msec being the primary's and secondary's masseso Dividing fll29p by the primary's 
principal moment of inertia ^ Mprimary , we obtain the contribution that this component of 
the torque brings into the angular deceleration rate 6 : 



9 = JC 
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ye^ + ^e^ + 0(e^) 
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(129) 
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+ O(zVQ) + 0(Q-3) + 0(ae2g-V/^) (130a) 
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the factor /C being given by 



(130b) 
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where ^ is a multiplier emerging in the expression ^ Mprimary R^ for the primary's principal 
moment of inertia (^ = 2/5 for a homogeneous sphere). 



For 9 :s> 2n, all the modes ^^2209 negative, so w^^^^ 
we process (|127p . for q — 1, like 



X220q ■ Then, keeping in mind that njQ <^ 1 . 



(131) 
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and similarly for other 9 7^ 0, and then plug the results into (|124p . This renders us (|128p . 
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Impq 



To arrive at the right-hand side of (|129p . we recalled that Ximpq ^^impq 
O(e^) = |e/™„„l + OiQ^^), according to formula 



Impq I 



l + ^-+Oiin/9r) 



and that Q, ^ 

^ Impq 
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In the special case of a = — 1 , the above expressions enjoy agreement with the appropriate 
result stemming from the corrected MacDonald model - except that our (11281) and (11301) contain 
At22oo ' '^2200 ' X2200 instead of At , Q , and x standing in formulae (44 - 47) from Williams 
& Efroimsky (2012). Formula ( I130bp tells us that the secular part of the tidal torque vanishes 



for 

0-n = -6ne^a , (132) 

which coincides, for a = —1 , with the result obtained in Rodriguez et al. (2008, eqn. 2.4) and 
Williams & Efroimsky (2012, eqn. 49). This coincidence however should not be taken at its face 
value, because it is occasional or, possibly better to say, exceptional. 

Formulae (11281 - 11291) were obtained by insertion of the expressions for the eccentricity func- 
tions and the phase lags into (I124p . and by assuming that n ^ \6\ . The latter caveat is a crucial 
element, not to be overlooked by the users of formulae (11281 - 11291) and of their corollary (11301 - 
I13ip for the tidal deceleration rate. 

The case of a = —1 is special, in that it permits derivation of (11281 - [TS^ without assuming 

that n \ \ . However for a > —1 the condition n ^ | | remains mandatory, so formulae 

(11281 - 11311) become inapplicable when reduces to values of about several n . 

Although formulae (I128ap and (I130ap contain an absolute error 0{ae'^Q^^n/6) , this does 

not mean that for comparable to n the absolute error becomes 0{ae'^Q~^) and the relative 

one becomes 0{ae^) . In reality, for comparable to n, the entire approximation falls apart, 
because formulae (11261 - I127P were derived from expression (11250 , which is valid for Q ^ 1 only 
(unless a = — 1 ). So these formulae become inapplicable in the vicinity of a commensurability. 
By ignoring this limitation, one can easily encounter unphysical paradoxesl^ 

Thence, in all situations, except for the unrealistic rheology a = —1 , limitations of the 
approximation (11281 - 11311) should be kept in mind. This approximation remains acceptable for 
n \6\ , but becomes misleading on approach to the physically- interesting resonances. 

9.2 An oversight in a popular formula 

The form in which our approximation (11281 - 11311) is cast may appear awkward. The formula for 

the despinning rate 6 is written as a function of 6 and n , multiplied by the overall factor /C . 

This form would be reasonable, were /C a constant. That this is not the case can be seen from 

the presence of the multiplier -tt;^ — = — on the right-hand side of ( 11311) . 

^2200 2\e-n\ 

Still, when written in this form, our result is easy to juxtapose with an appropriate formula 
from Correia & Laskar (2004, 2009). There, the expression for the despinning rate looks similar 
to ours, up to an important detail: the overall factor is a constant, because it lacks the said 
multiplier ^r;^ — . The multiplier was lost in those two papers, because the quality factor was 

A2200 

introduced there as l/(nAt), see the line after formula (9) in Correia & Laskar (2009). In 
reality, the quality factor Q should, of course, be a function of the tidal frequency x 5 because 
Q serves the purpose of describing the tidal damping at this frequency. Had the quality factor 
been taken as l/(x At) , it would render the corrected MacDonald model (a = — 1 ), and the 
missing multiplied would be there. Being unphysical 1^ the model is mathematically convenient, 
because it enables one to write down the secular part of the torque as one expression, avoiding 



For example, in the case of a > — 1 , formulae (|125l - [TTl\ render infinite values for Aij^^^^ and Cj^p^ on 
crossing the commensurability, i.e., when Wj^^^ goes through zero. 

To be exact, the model is unphysical everywhere except in the closest vicinity of the resonance - see formulae 
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the expansion into a scries (Williams & Efroimsky 2012). The model was pioneered by Singer 
(1968) and employed by Mignard (1979, 1980), Hut (1981) and other authors. 

Interestingly, in the special case of the 3:2 spin-orbit resonance, we have x = n . Still, 
the difference between x ^"^^ ^ the vicinity of the resonance may alter the probability of 
entrapment of Mercury into this rotation mode. The difference between x ^-iid n becomes even 
more considerable near the other resonances of interest. So the probabilities of entrapment into 
those resonances must be recalculated. 

10 Conclusions 

The goal of this paper was to lay the ground for a reliable model of tidal entrapment into 
spin-orbital resonances. To this end, we approached the tidal theory from the first principles 
of solid-state mechanics. Starting from the expression for the material's compliance in the time 
domain, we derived the frequency- dependence of the Fourier components of the tidal torque. The 
other torque, one caused by the triaxiality of the rotator, is not a part of this study and will be 
addressed elsewhere. 

• We base our work on the Andrade rheological model, providing arguments in favour of its 
applicability to the Earth's mantle, and therefore, very likely, to other terrestrial planets 
and moons. The model is also likely to apply to the icy moons (Castillo-Rogez et al. 
2011). We have reformulated the model in terms of a characteristic anelastic timescale 

(the Andrade time). The ratio of the Andrade time to the viscoelastic Maxwell time, 
= t^/Tj^j, serves as a dimensionlcss free parameter of the rheological model. 

The parameters r^, r^j,^, ( cannot be regarded constant, though their values may be changing 
very slowly over vast bands of frequency. The shapes of these frequency-dependencies may 
depend on the dominating dissipation mechanisms and, thereby, on the magnitude of the 
load, as different damping mechanisms get activated under different loads. 

The main question here is whether, in the low-frequency limit, anelasticity becomes much 
weaker than viscosity. (That would imply an increase of and C as the tidal frequency 
X goes down.) The study of ices under weak loads, with friction caused mainly by lattice 
diffusion (Castillo-Rogez et al. 2011, Castillo-Rogez & Choukroun 2010) has not shown 
such a decline of anelasticity. However, Karato & Spetzler (1990) point out that it should 
be happening in the Earth's mantle, where the loads are much higher and damping is caused 
mainly by unpinning of dislocations. According to Ibid., in the Earth, the decrease of the 
role of anelasticity happens abruptly as the frequency falls below the threshold Xq yr~^ • 
We then may expect a similar switch in the other terrestrial planets and the Moon, though 
there the threshold may be different as it is sensitive to the temperature of the mantle. 
The question, though, remains if this statement is valid also for the small bodies, in which 
the tidal stresses are weaker and dissipation is dominated by lattice diffusion. 

• By direct calculation, we have derived the frequency dependencies of the factors sin 
emerging in the tidal theory. Naturally, the obtained dependencies of these factors upon 
the tidal frequency Ximpq (°^' P^^ more exactly, upon the tidal mode w^^^p^ ) mimic 
the frequency-dependence of the imaginary part of the complex compliance. They scale as 
~ X~" with < q; < 1 , at higher frequencies; and look as ~ lower frequencies. 
However in the zero-frequency limit the factors A;^ sin demonstrate a behavior inherent 
in the tidal lagging and absent in the case of lagging in a sample: in a close vicinity of the 
zero frequency, these factors (and the appropriate components of the tidal torque) become 
linear in the frequency. This way, sin first reaches a finite maximum, then decreases 
continuously to nil as the frequency approaches to zero, and then changes its sign. So 
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the resonances are crossed continuously, with neither the tidal torque nor the tidal force 
diverging there. For example, the leading term of the torque vanishes at the synchronous 
orbit. 

This continuous traversing of resonances was pointed out in a heuristic manner by Efroimsky 
& Williams (2009). Now we have derived this result directly from the expression for the 
compliance of the material of the rotating body. Our derivation, however, has a problem 
in it: the frequency, below which the factors sine^ and the appropriate components 
of the torque change their frequency-dependence to linear, is implausibly low (lower than 
10~^°Hz, if we take our formulae literally). The reason for this mishap is that in our 
formulae we kept using the known value of the Maxwell time r^,,^ all the way down to the 
zero frequency. Possible changes of the viscosity and, accordingly, of the Maxwell time in 
the zero-frequency limit may broaden this region of linear dependence. 

• We have offered an explanation of the unexpected frequency-dependence of dissipation 
in the Moon, discovered by LLR. The main point of our explanation is that the LLR 
measures the tidal dissipation whose frequency-dependence is different from that of the 
seismic dissipation. Specifically, the "wrong" sign of the exponent in the power dissipation 
law may indicate that the frequencies at which tidal friction was observed were below the 
frequency at which the lunar k2 sin 62 has its peak. Taken the relatively high frequencies 
of observation (corresponding to periods of order month to year), this explanation can be 
accepted only if the lunar mantle has a low mean viscosity. This may be the case, taken 
the presumably high concentration of the partial melt in the low mantle. 

• We have developed a detailed formalism for the tidal torque, and have singled out its 
oscillating component. 

The studies of entrapment into spin-orbit resonances, performed in the past, took into 
account neither the afore-mentioned complicated frequency- dependence of the torque in 
the vicinity of a resonance, nor the oscillating part of the torque. We have written down 
a concise and practical expression for the oscillating part, and have raised the question 
whether it may play a role in the entrapment and libration dynamics. 
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Appendix. 



A Continuum mechanics. 

A celestial-mechanician's survival kit 

Appendix |X] offers an extremely short introduction into continuum mechanics. The standard 
material, which normally occupies large chapters in books, is compressed into several pages. 

Subsection lA.ll presents the necessary terms and definitions. Subsections IA.2I explains the 
basic concepts employed in the theory of stationary deformation, while subsection IA.3I explains 
extension of these methods to creep. These subsections also demonstrate the great benefits stem- 
ming from the isotropy and incompressibility assumptions. Subsection IA.4I introduces viscosity, 
while subsection IA.5I offers an example of how elasticity, viscosity get combined with hereditary 
reaction, into one expression. Subsection IA.6I renders several simple examples. 

A.l Glossary 

We start out with a brief guide elucidating the main terms employed in continuum mechanics. 

• Rheology is the science of deformation and flow. 

• Elasticity: This is the most trivial reaction - instantaneous, linear, and fully reversible 
after stressing is turned off. 

• Anelasticity: While still linear, this kind of deformation is not necessarily instantaneous, 
and can demonstrate "memory", both under loading and when the load is removed. Im- 
portantly, the term anelasticity always implies reversibility: though with delay, the original 
shape is restored. Thus an anelastic material can assume two equilibrium states: one is the 
unstressed state, the other being the long-term relaxed state. Anelasticity is characterised 
by the difference in strain between the two states. It is also characterised by a relaxation 
time between these states, and by its inverse - the frequency at which relaxation is most ef- 
fective. The Hohenemser-Prager model, also called SAS (Standard Anelastic Sohd), renders 
an example of anelastic behaviour. 

Anelasticity is an example of but not synonym to hereditary reaction. The latter includes 
also those kinds of delayed deformation, which are irreversible. 

• Inelasticity: This term embraces all kinds of irreversible deformation, i.e., deformation 
which stays, fully or partially, after the load is removed. 

• Unelasticity ( = N anelasticity): These terms are very broad, in that they embrace any be- 
haviour which is not elastic. In the literature, these terms are employed both for recoverable 
(anelastic) and unrecoverable (inelastic) deformations. 

• Plasticity: Some authors simply state that plastic deformation is deformation which is 
irreversible - a very broad and therefore useless definition which makes plasticity sound 
like a synonym to inelasticity. 

More precisely, plastic is a stepwise behaviour: no strain is observed until the stress a 
reaches a threshold value (called yield strength), while a steady flow begins as the 
stress reaches the said threshold. Plasticity can be either perfect (when deformation is 
going on without any increase in load) or with hardening (when increasingly higher stresses 
are needed to sustain the flow). It is always irreversible. 
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In real life, plasticity shows itself in combination with elasticity or/and viscosity. Models 
describing these types of behaviour are called elastoplastic, viscoplastic, and viscoelasto- 
plastic. They are all inelastic, in that they describe unrecoverable changes of shape. 

• Viscosity: Another example of inelastic, i.e., irreversible behaviour. A viscous stress is 
proportional to the time derivative of the viscous strain. 

• Viscoelasticity: The term is customarily applied to all situations where both viscous and 
elastic (but not plastic) reactions are observed. One may then expect that the equations 
interrelating the viscoelastic stress to the strain would contain only the viscosity coefficients 
and the elastic moduli. However this is not necessarily true, as some other empirical 
constants may show up. For example, the Andrade model (IHTj) contains an elastic term, 
a viscous term, and an extra term responsible for a hereditary reaction (the "memory"). 
Despite the presence of that extra term, the Andrade creep is still regarded viscoelastic. So 
it should be understood that viscoelasticity is, generally, more than just viscosity combined 
with elasticity. One might christen such deformations "viscoelastohereditary" , but such a 
term would sound awkward 

On many occasions, complex viscoelastic models can be illustrated with an infinite set of 
viscous and elastic elements. These serve to interpret the hereditary terms as manifestations 
of viscosity and elasticity only. While illustrative, these schemes with dashpots and springs 
have their limitations and should not be taken too literally. In some (not all) situations, 
the hereditary terms may be interpreted as time- dependent additions to the steady-state 
viscosity coefficient, the Andrade model being an example of such situation. 

• Viscoplasticity: These are all models wherein both viscosity and plasticity are present 
in some combination. In these situations, higher stresses have to be applied to increase 
the deformation rate. Just as in the case of viscoelasticity, viscoplastic models may, in 
principle, incorporate extra terms standing for hereditary reaction. 

• Elastoviscoplasticity ( = Viscoelastoplasticity) : The same as above, though with elasticity 
present. 

• Hereditary reaction: While the term is self-explanatory, it would be good to limit its use to 
effects other than viscosity. In expression fH6|l for the stress through strain, the distinction 
between the viscous and hereditary reactions is clear: while the viscous part of the stress is 
rendered instantaneously by the delta-function term of the kernel, the hereditary reaction 
is obtained through integration. In expression ( HOj) for the strain through stress, though, 
the viscous part shows up, under the integral, in sum with the other delayed terms - see, 
for example, the Andrade model (IHTIl . This is one reason for which we shall use the term 
hereditary reaction in application to delayed behaviour different from the pure viscosity. 
Another reason is that viscous flow is always irreversible, while a hereditary reaction may 
be either irreversible (inelastic) or reversible (anelastic). 

• Creep: Widely accepted is the convention that this term signifles a slow-rate deformation 
under loads below the yield strength . 

Numerous authors, though, use the oxymoron plastic creep, thereby extending the applica- 
bility realm of the word creep to any slow deformation. 

Here we shall understand creep in the former sense, i.e., with no plasticity involved. 

' Sometimes the term elastohereditary is used, but not viscoelastohereditary or elastoviscohereditary. 
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It would be important to distinguish between viscoelastic deformations, on the one hand, and 
viscoplastic (or, properly speaking, viscoelastoplastic) deformations on the other hand. Plasticity 
shows itself at larger stresses and is, typically, nonlinear. It comes into play when the linearity 
assertion fails. For most minerals, this happens when the strain approaches the threshold of 
10~^ . Although it is possible that this threshold is transcended in some satellites (for example, 
in the deeper layers of the Moon), we do not address plasticity in this paper. 



A. 2 Stationary linear deformation of isotropic incompressible media 

In the linear approximation, the tensor of elastic stress, § , is proportional to the differences in 
displacement of the neighbouring elements of the medium. These differences are components of 
the tensor gradient V ® u , where u is the displacement vector. 

The tensor gradient can be decomposed, in an invariant way, into its symmetric and antisym- 
metric parts: 



V (g) u 



1 
2 



(V®u) + (V®u) 



1 
2 



(V 



u 



(V®u) 



(133) 



The decomposition being invariant, the two parts should contribute into the stress independently, 
at least in the linear approximation. However, as well known (Landau & Lifshitz 1986), the 
antisymmetric part of (11331) describes the displacement of the medium as a whole and thus 



brings nothing into the stress. This is why the linear dependence is normally written as 



§ = B U 



(134) 



where B is a four-dimensional tensor having 3'^ = 81 components, while the strain tensor 



U 



u 



+ (V 



u 



(135) 



is the symmetric part of the tensor gradient. Its components are related to the displacement 



vector u through u^^ ^ 1 (^|^ + 



Although the matrix B is comprised of 81 empirical constants, in isotropic materials the 
description can be reduced to two constants only. To see this, recall that the expansion of 



the strain into a part with trace and a traceless part, U = ^ISpU + ^ 



U 



ISpU 



IS 



invariant. Here the trace of U is denoted with SpU 



Ur 



summation over repeated indices 



being implied. The notation I stands for the unity matrix consisting of elements 5.^^ ■ 

In an isotropic medium, the elastic stress must be invariantly expandable into parts propor- 
tional to the afore-mentioned parts of the strain. The first part of the stress is proportional, with 
an empirical coefficient 3 K ,to the trace part ^ I Sp U of the strain. The second part of the stress 

will be proportional, with an empirical coefficient 2 /i , to the traceless part ^U — -^ISpUj of 
the strain: 



K I SpU + 2i2{V 



ilSpU 



j9l + 2/i( U - ^ISpU 



(136a) 



or, in Cartesian coordinates: 



p 5^1, + 2 n \ u. 



2 ^yv 'U'aa 



(136b) 
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where 

p = - ^ Sp§ = -KSpV (137) 
is the hydrostatic pressure. Thus the elastic stress gets decomposed, in an invariant way, as: 

(e) _ (e) (e) 

^ ~ ^volumetric ^ deviatoric ' (138) 

where the volumetric elastic stress is given by 



(£) „ „ „ „ 1 

3 



.olumetr^c - K I SpV = I - Sp S = - pi , (139) 



while the deviatoric elastic stress is: 

sLato„c = 2/x(u- ilSpU^ . (140) 
Inverse to ([I36a| - [HeH) are the following expressions for the strain tensor: 

U=— ISp§+— f§--ISp§J (141a) 

and 

y A Z fi \ 6 J 

1 (e) 

where the term with trace, g-^ I Sp S , is called the volumetric strain, while the traceless 

1 /(e) 1 (e)\ 

term, ^-^ (§— ^ISp§j,is named the deviatoric strain. The quantity K is called the ^ti/A; 

modulus, while is called the shear modulus. 

Expressions f ll36p trivially entail the following interrelation between traces: 

(e) (e) 

Sp §= 3/^SpU or, in terms of components: craa= 3 K Uaa ■ (142) 

As demonstrated in many books (e.g., in Landau & Lifshitz 1986), the trace of the strain is 
equal to the relative variation of the volume, experienced by the material subject to deformation: 
Uaa = V • u = , where u is the displacement vector. In the no-compressibility 

approximation, the trace of the strain and, according to fll42p . that of the stress become zero. 
Then, in the said approximation, the hydrostatic pressure (11371) and the volumetric elastic stress 
(11391) become nil, and all we are left with is the deviatoric elastic stress (and, accordingly, the 
deviatoric part of the strain). Formulae f ll36p and fll4ip get simplified to 

(e) 

§ = 2/iILJ , which is the same as a^^ = , (143) 



and to 



(e) 

2 U = J § , which is the same as 2u^^ = J a , (144) 



the quantity J = 1/ ^ being called the compliance of the material. 



41 



A. 3 Evolving linear deformation of isotropic incompressible isotropic 
media. Hereditary reaction 

Equations (11341 - I144p were written for static deformation, so each of these equations can be 
assumed to connect the strain and the elastic stress taken at the same instant of time (for a 
static deformation their values stay constant anyway). 

Extension of this machinery is needed when one wants to describe evolving deformation of 
materials with "memory". Thence the four-dimensional tensor B becomes a linear operator B 
acting on the strain tensor function as a whole. To render the value of the stress at time t , the 
operator will "consume" as arguments all the values of strain over the interval t' G ( — oo , t ] : 

(h) 

§ [t) = B(t) U . (145) 

Thus B will be an integral operator, with integration going from t' = — oo through t' = t . 

In the static case, the linearity guaranteed elasticity, i.e., the abihty of the body to regain 
its shape after the loading is turned off: no stress yields no strain. In a more general situation 

of materials with "memory" , this ability is no longer retained, as the material may demonstrate 

(h) 

creep. This is why, in this section, the stress is called hereditary and is denoted with § . 

Just as in the stationary case, we wish the properties of the medium to remain isotropic. 
As the decomposition of the strain into the trace and traceless parts remains invariant at each 
moment of time, these two parts will, separately, generate the trace and traceless parts of the 
hereditary stress in an isotropic medium. This means that, in such media, the four- dimensional 
tensor operator B gets reduced to two scalar linear operators K and /i : 

(h) {h) (h) _ / 1 \ 

S = ^ .olumetr^c + ^ de.^ator^c = K I SpV + 2 i^V - -I SpV j , (146) 

where both K and jl are integral operators acting on the tensor function u^u{t') as a whole, 
i.e., with integration going from t' = — oo through t' = t. 

If we also assume that, under evolving load, the medium stays incompressible, the trace of the 
strain, Uaa , will stay zero. An operator generalisation of fll42p now reads: (Taait) = 3 K{t) Uaa ■ 
Under a reasonable assumption of being nil in the distant past, this integral operator renders 
o'aa = at all times. This way, in a medium that is both isotropic and incompressible, we have: 

(h) {h) 

U = ^deviatoric ^ud, accordingly: § = . (147) 

Then the time-dependent analogues to formulae ( 1143P and ( 1144p will be: 

(h) 

§ (t) = 2fi{t)V (148) 

and 

2U(t) = J{t) § , (149) 

where the compliance J, too, has been promoted to operatorship and crowned with a caret. 

Formula fll48p tells us that in a medium, which is both isotropic and incompressible, relation 
between the stress and strain tensors can be described with one scalar integral operator fi only, 
the complementary operator J being its inverse. (Here the adjective "scalar" does not imply 
multiplication with a scalar. It means that the operator preserves its functional form under a 
change of coordinates.) 

Below we shall bring into the picture also the viscous component of the stress, a compo- 
nent related to the strain through a four-dimensional tensor whose 3^ = 81 components are 
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differential operators. In that case too, the isotropy of the medium will enable us to reduce the 
81-component tensor operator to two differential operators transforming as scalars. Besides, the 
incompressibility of the medium makes the viscous stress traceless. Thus it will turn out that, in 
an isotropic and incompressible medium, the viscous component of the stress can be described 
by only one scalar differential operator - much like the elastic and hereditary parts of the stress. 
(Once again, "scalar" means: indifferent to coordinate transformations.) 

Eventually, the elastic, hereditary, and viscous deformations will be united under the auspices 
of a general viscoelastic formalism. In an isotropic medium, this combined formalism will be 
reduced to two integro-differential operators only. In a medium which is both isotropic and 
incompressible, the formalism will be reduced to only one scalar integro-differential operator. 



A. 4 The viscous stress 



(e) (v) 

While the elastic stress § is linear in the strain, the viscous stress § is linear in the first 
derivatives of the components of the velocity with respect to the coordinates: 



(v) 



(V 



(150) 



where A is the so-called viscosity tensor, V ® v is the tensor gradient of the velocity. The 
velocity of a fluid parcel relative to its average position is connected to the displacement vector 
u through V = du/dt . 

The tensor gradient of the velocity can be expanded, in an invariant way, into its antisym- 
metric and symmetric parts: 



V®v = f2 + E , 
where the antisymmetric part is furnished by the vorticity tensor 



(V 



(V 



while the symmetric part is given by the rate-of-shear tensor 



E 



1 
2 



(V®v) + (V®v) 



The latter is obviously related to the strain tensor through 



(151) 



(152) 



(153) 



(154) 



It can be demonstrated (e.g.. Landau & Lifshitz 1987) that the antisymmetric vorticity tensor 
describes the rotation of the medium as a wholj^ and therefore contributes nothing to the 
stress@ For this reason, the viscous stress can be written as 



(v) 



.E 



d_ 

di 



U 



(155) 



The matrix A is four-dimensional and contains 3'^ = 81 components. Just as the matrix B 
emerging in expression (11341) for the elastic stress, the matrix A can be reduced, in an isotropic 



^21 



This is why this tensor's components coincide with those of the angular velocity u5 of the body. For example, 
= ^) '^oi^'^ides with ^3 . 

^^Since expansion (I15ip of the tensor gradient into the vorticity and rate-of-shear tensors is invariant, then so 
is the conclusion about the irrelevance of the vorticity tensor for the stress picture. 
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medium, to only two empirical constants. To see this, mind that the rate-of-shear tensor can be 
decomposed, in an invariant manner, into two parts: 



E=^IV-v + (e-^IV-v) . (156) 



where the rate- of- expansion tensor 



^ I V ■ V (157) 



is diagonal and has a trace, while the combination 



E- ilV.v= i 



(V®v) + (V®v)^ -ilV-v (158) 

o 



is symmetric and traceless. These two parts contribute linearly proportional inputs into the 
stress. The first input is proportional, with an empirical coefficient 3 C , to the rate-of-expansion 
term, while the second input into the stress is proportional, with an empirical coefficient 2 77 , 
to the symmetric traceless combination: 

§= 3C^IV-v + 2r/ (^E - ilV-v^ = (ISpU) + 2r/|^ (^U - ^ I SpU^ . (159) 

Here we recalled that V-v = ^V-u = ^ Uaa = ^ ^ • Since Sp U is equal to the 

volume variation Jy^^ experienced by the material, we see that the ffist term in fll59p is 
volumetric, the second being deviatoric. 

The quantity 77 is called the first viscosity or the shear viscosity. The quantity ( is named 
the second viscosity or the bulk viscosity. 



A. 5 An example of approach to viscoelastic behaviour 

In this subsection, we shall consider one possible approach to description of viscoelastic regimes. 
As we mentioned in subsection lA.ll the term viscoelasticity covers not only combinations of elas- 
ticity and viscosity, but can also include other forms of delayed reaction. So the term viscoelastic 
is customarily used as a substitution for too long a term viscoelastohereditary. 

One possible approach would be to assume that the elastic, hereditary, and viscous stresses 
simply sum up, and that each of them is related the the same strain U : an 

{total) (e) [h) [v) f ^ Q \ 

+ B + A— U , (160a) 



or simply 



dt 



{total) „ ^ ^ Q 

§ = BU , where B = B + B + A— , (160b) 

at 



8 

where the three operators - the integral operator B , the differential operator A ^ , and the 
operator of multiplication by matrix B - comprise an integro-differential operator B . 

In an isotropic medium, each of the three matrices, B , A ^ , and B , includes two terms only. 
This happens because in such a medium each of the three parts of the stress gets decomposed 
invariantly into its deviatoric and volumetric components: The elastic stress becomes: 



S = ^tolu.^etr^c + ,e.^ator^c =3irQlSpu)+2^(u-ilSpu) , (161) 
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with K and /i being the bulk elastic modulus and the shear elastic modulus, correspondingly, I 
standing for the unity matrix, and Sp denoting the trace of a matrix: Sp U = ^ - Uu . 
The hereditary stress becomes: 



(h) (h) (h) ~ / 1 \ / I 

-ISpUj + 2/i f U- - 



.olume,r^c + ^ de.^ator^c = 3 ( - I Sp U ) + 2 /i ( U - - I Sp U ) , (162) 



where K and are the hulk-modulus operator and the shear-modulus operator, accordingly. 
The viscous stress acquires the form: 



S = S.oZ.„.e*„c + Sde«a*o„c= ^ C ^ (^3 I Sp U j + 2 7/ - U - - I Sp U j , (163) 

the quantities C, and r] being termed as the hulk viscosity and the shear viscosity, correspondingly 
The term ^ I Sp U is called the volumetric part of the strain, while U — ^ I Sp U is called 

the deviatoric part. Accordingly, in expressions (11611 - [T63|) for the stresses, the pure-trace terms 

are called volumetric, the other term being named deviatoric. 

The total stress, too, can now be split into the total volumetric and the total deviatoric parts: 

(e) (v) (h) 

§ § § 

^ 



(total) /(e) (e) \ / (v) (v) \ f {h) (h) 



volumetric ' ^ deviatoric ; I volumetric ' ^ deviatoric ; I volumetric ' ^ deviatoric 



volumetric deviatoric 



■(e) (v) (h) \ fie) (v) (h) 

volumetric ' volumetric ' volumetric j \ deviatoric ' deviatoric ' deviatoric 



3K + 3k + SC^] (^ISpu") + (2^1 + 2/i + 2r]^^ (v - ^ISpU ) (164a) 



= SxQlSpU^ + 2 /i - ^ISpU^ , (164b) 

where 

d d 
K = K + K-^C— and fl = jj, + fM + ri— . (165) 
ot at 

As expected, a total linear deformation of an isotropic material can be described with two integro- 
differential operators, one acting on the volumetric strain, another on the deviatoric strain. 

If an isotropic medium is also incompressible, the relative change of the volume vanishes: 
Sp U = . Accordingly, the volumetric part of the strain becomes nil, and so do the volumetric 
parts of the elastic, hereditary, and viscous stresses. For such media, we end up with a simple 
relation which includes only deviators: 

[total) ^ ^(e) (h) (v) _ d 

S ~ ^deviatoric ~ '^deviatoric + '^deviatoric + '^ deviatoric ~ 2/iU-|-2/iU-|-2?7— U (166) 



or simply: 

(total) 



deviatoric 



2/iU , (167) 
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where U contains only a deviatoric part, while 



d 

/i = /x + /i + ?7— (168) 

is the total, integro-differential operator, which is mapping the preceding history and the present 
rate of change of the strain to the present value of the stress. 

It should be reiterated that the above approach is based on the assertion that the elastic, 
viscous, and hereditary stresses sum up, and that all three are related to the same total strain. A 
simple example of this approach, called the Kelvin- Voigt model, is rendered below in subsection 

A different approach would be to assume that the strain consists of three distinct parts - 
elastic, hereditary, and viscous - and that these components are related to the same overall 
stress. A simple example of this treatment, termed the Maxwell model, is set out in subsection 
IA.6.41 A more complex example of this approach is furnished by the Andrade model presented 
in subsection 15.31 Another way of combining elasticity and viscosity (with no hereditary reaction 
involved) is implemented by the Hohenemser-Prager (SAS) model explained in subsection IA.6.51 
below. 

A. 6 Examples of viscoelastic behaviour with no hereditary reaction 
A. 6.1 Elastic deformation 

The truly simplest example of deformation is elastic: 

(e) (e) ^ ^ 

§=2/iU , U = J§ , (169) 
where fj, and J are the unrelaxed rigidity and compliance: 

^ = ^(0) , J = J(0) , /i = 1/J . (170) 
In the frequency domain, this relation assumes the same form as it would in the time domain: 

^^Ax) = '^l^u^M , 2m^^(x) = J^^u{x) ■ (171) 
A. 6. 2 Viscous deformation 

The next example is that of a purely viscous behaviour: 

i^) d , , 

§= 2r/^U . (172) 

It is straightforward from fll72p and fl271) that, in this regime, the Fourier components of the 
stresQ are connected to those of the strain through 

o^iv{.x) = Hx) u^uix) , 2m^^(x) = Jix)^'yuix) , (173) 
where the complex rigidity and the complex compliance are given by 

- i 

IX = irix , J = ■ (174) 

VX 



Although we no longer spell it out, the word stress everywhere means: deviatoric stress, as we agreed to 
consider the medium incompressible. 
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A. 6. 3 Viscoelastic deformation: a Kelvin- Voigt material 



The Kelvin- Voigt model, also called the Voigt model, can be represented with a purely viscous 
damper and a purely elastic spring connected in parallel. Subject to the same elongation, these 
elements have their forces summed up. This illustrates the situation where the total, viscoelastic 
stress consists of a purely viscous and a purely elastic inputs called into being by the same strain: 



(lie) (v) (e) {v) (e) 

§ = § + § , while U = U = U 



Then the total stress reads: 



§ = ( 2/i + 2r]^ 1 U 



which is often presented in the form of 
with the so-called Voigt time defined as 



(175) 



(176a) 



;i76b) 



(177) 



Comparing (11761) with fH6l) . we understand that the kernel of the rigidity operator for the Kelvin- 
Voigt model can be written down as 



fi(t — t') = jl + T] 6(t — t' 

Suppose the strain is varying in time as 
so that 



2/i 



1 — exp 



t - to 



e(t - to) 



^71^ (^) = t; exp 

2/i 



t - to 



e(t - to) 



(178) 



(179) 



(180) 



Then insertion of f ll78p and fll80p into (??) or, equivalently, insertion of fll79p into f ll76ap 
demonstrates that this strain results from a strescll 

a^^it) = (To e(t - to) • (181) 

It would however be a mistake to deduce from this that the compliance function is equal to 



1 — exp 



t - t' 



, even though such a misstatement is sometimes made in the lit- 
erature. This expression furnishes the compliance function only in the special situation of a 
Heaviside-step stress fll8ip . but not in the general case. 



For example, plugging of ()178p and ()180p into (??) leads to: 



a{t) 



t '= t 



[fi + 1] 6{t — t')] — — exp 



t' - to 



Q{t' - to)dt' 



Co 





for t > to 
for t < to , 



which is simply cto ©(t — to) . 
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As can be easily shown from ( l27j) . in the frequency domain model f ll76p reads as (11731) . except 
that the complex rigidity and the complex compliance are now given by 

fi = + zxr^) , J = • (182) 

Recall that, for brevity, here and everywhere we write /i and J instead of yu(0) and J(0) . 

The Kelvin- Voigt material becomes elastic in the low- frequency limit, and viscous in the 
high-frequency limit. 



A. 6.4 Viscoelastic deformation: a Maxwell material 

The Maxwell model can be represented with a viscous damper and an elastic spring connected 
in series. Experiencing the same force, these elements have their elongations summed up. This 
example illustrates the situation where the total, viscoelastic strain consists of a purely viscous 
and a purely elastic contributions generated by the same stress § : 

{v) (e) (e) (e) (v) Q (v) 

U = U + U , where §=2/iU and §=2r/— U • (183) 

at 

Since in the Maxwell regime both contributions to the strain are generated by the same stress 

(ve) (v) (e) 



formula (11831) can be written down as 



184) 



1LJ= 7^ S +7^ § (185a) 



or, in a more conventional form: 



+ — § = 2/i U , (185b) 



with the so-called Maxwell time introduced as 

r„ = . (186) 

Although formally the Maxwell time is given by an expression mimicking the definition of the 
Voigt time, the meaning of these times is different. 

Comparing (I185p with the general expression for the compliance operator, we see that, 
for the Maxwell model, the compliance operator in the time domain assumes the form: 



Jit - t') 



J +{t - t')- 
V 



e{t - t') , (187) 



where J = 1/yU . In the frequency domain, (I185p can be written down as (I173p . with the complex 
rigidity and compliance given by 

Kx) = /i Y^T^ ' ^"(^) = ^ f 1 - = J - — ■ (188) 

1 + ^xm V x^Ai / XV 

Clearly, such a body becomes elastic in the high-frequency limit, and becomes viscous at low 
frequencies (the latter circumstance making the Maxwell model attractive to seismologists). 
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A. 6. 5 Viscoelastic deformation: the Hohenemser-Prager (SAS) model 



An attempt to combine tlie Kelvin- Voigt and Maxwell models leads to the Hohenemser-Prager 
model, also known as the Standard Anelastic Solid (SAS): 

r,,§ + § = 2/i (^U + vil) , (189) 

In the limit of r^^ — )■ , this model approaches the one of Kelvin and Voigt (and acquires 
the meaning of the Voigt time). 

A transition from the SAS to Maxwell model, however, can be achieved only through re- 
definition of parameters. One should set: 2 yU — )■ and — )■ oo , along with 2/ir^,, — t- 21] . 
Then f ll89p will become fllSSp . with r^^ playing the role of the Maxwell time. 

In the frequency domain, fll89p can be put into the form of fll73p . the complex rigidity and 
the complex compliance being expressed through the parameters as 

p = ,'-±^ . J=J\^ (190) 



[^v-^m)X 

is related to its maximal value through 



This entails: tan 5 = Zm[fj]/7le[^] = y~t ^ -2 ? whence it is easy to show that the tangent 



tan 5 = 2 [tan 5]^^^ ^ + ' ^^'^'''^ ^ ^ ^^^^ ' 

This is the so-called Debye peak, which is indeed observed in some materials. 

To prove that the SAS solid is indeed anelastic, one has to make sure that a Heaviside step 
stress Q{t') entails a strain proportional to l—exp(— Ft) , and to demonstrate that a predeformed 
sample subject to stress 9(— t') regains it shape as exp(— Ft) , where the relaxation constant F 
is positive. In subsection IC.3I we shall do this for a SAS sphere. 



B Interconnection between the quality factor and the phase 
lag 

The power P exerted by a tide-raising secondary on its primary can be written as 

p = - [ p V ■ VWd^x (191) 



p , V , and W signifying the density, velocity, and tidal potential in the small volume (Px 

of the primary. The mass-conservation law V ■ (p V) + ^ = enables one to shape the 
dot-product into the form of 

pV-VW = V-{pVW)-pWV-V-VWVp . (192) 

— * 

Under the realistic assumption of the primary's incompressibility, the term with V ■ V may be 
omitted. To get rid of the term with Vp, one has to accept a much stronger approximation of 
the primary being homogeneous. Then the power will be rendered by 

P = - j V ■ {p V W)£x = - j pW V ■ a dS , (193) 

n being the outward normal and dS being an element of the surface area of the primary. This 
expression for the power (pioneered, probably, by Goldreich 1963) enables one to calculate the 
work through radial displacements only, in neglect of horizontal motion. 
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^ = = 4 • (194) 



Denoting the radial elevation with , we can write the power per unit mass, V = P/M , a.s: 

V n - (- 

9r / \ dr J dt 

A harmonic external potential 

W = Wo cos{u,^^^t) , (195) 
applied at a point of the primary's surface, will elevate this point by 

Wo Wo 

C = h2— cos^^p^t - e^^p^) = h2— cos^^p^t - w^^p^At^^p^) , (196) 

with g being the surface gravity acceleration, and /i2 denoting the Love number. 

In formula fll96p . oji^p^ is one of the modes fll05p showing up in the Darwin-Kaula expansion 
dinSD- The quantity e^^^^ = w^^p, At;^pg is the corresponding phase lag, while At^^^^ is 
the positively defined time lag at this mode. Although the tidal modes 00^^^^ can assume any 
sign, both the potential W and elevation ( can be expressed via the positively defined forcing 
frequency Ximpq ~ I ^impq I absolute value of the phase lag: 

W = Wo cos(xt) , (197) 

C = h2— cos(xt - |e|) , (198) 
g 

subscripts Impq being dropped here and hereafter for brevity. 

The vertical velocity of the considered element of the primary's surface will be 

-7- = ~h2X~~ sin(xt — lei) = —h2X~— i^^T^Xi cos lei — cosx^ sin lei) . (199) 
dt g g 

Introducing the notation A = h2 ' write the power per unit mass as 

^ = f = cos(xt) sin(xt - |e|) , (200) 

and write the work w per unit mass, performed over a time interval (to , t) , as: 

w\^^= I V dt = A cos(xt) sin(xt — |e|)c/(xt) = A cos |e| / cosz sin2;(iz — A sin |e| / cos^ zdz 

= - j[ cos(2xt - |e|) + 2 X t sin |e| ] . (201) 

Being cyclic, the first term in (120 ip renders the elastic energy stored in the body. The second 
term, being linear in time, furnishes the energy damped. This clear interpretation of the two 
terms was offered by Stan Peale [2011, personal communication]. 

The work over a time period T = 27r/x is equal to the energy dissipated over the period: 

^Ulo =A^c,de = - ATI sin|e| . (202) 

It can be shown that the peak work is obtained over the time span from tt to |e| and assumes 
the value 



(work) A 

E ^ = - 

peak 2 



TT 

COS lei — sm |e| 



2 



(203) 
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whence the appropriate quahty factor is given by: 



tan 


e 




i-(5- 


e 


) 


1 tan 


e 



- AE , 

Q-i = = , . (204) 

^work (work) / ^ . . \ ^ ' 

2tx hi 

peak 

To calculate the peak energy stored in the body, we would note that the first term in (120 ip 
is maximal when taken over the span from xt = 7r/4 + |e|/2 through = 37r/4 + |e|/2 : 

E^T^^ = ^ , (205) 

peak 2 

and the corresponding quality factor is: 

- AE , 

= = sin|e| . (206) 

^energy {energy) i i ^ ' 

Zn h, 

peak 

Goldreich (1963) suggested to employ the span xt = (0, 7r/4) . The absolute value of the 
resulting power, denoted in Ihid. as E* , is equal to 

E* = - cos|e| (207) 

and is not the peak value of the energy stored nor of the work performed. Goldreich (1963) 
however employed it to define a quality factor, which we shall term Qcoidreich ■ '^^^^ factor is 
introduced via 



Q-\ = £^ = tan|e| . (208) 

^ Goldreich 2 Tf E* 



In our opinion, the quality factor Q^^ergy defined through f l206p is preferable, because the ex- 
pansion of tides contains terms proportional to kiiXimpq) ^^^^KXimpg) • Siucc the loug-established 
tradition suggests to substitute sine with l/Q , it is advisable to define the Q exactly as fl206p . 
and also to call it Qi , to distinguish it from the seismic quality factor (Efroimsky 2012). 



C Tidal response of a homogeneous viscoelastic sphere 
(Churkin 1998) 

This section presents some results from the unpublished preprint by Churkin (1998). We took 
the liberty of upgrading the notation^ and correcting some minor oversights. 



C.l A homogeneous Kelvin- Voigt spherical body 

In combination with the Correspondence Principle, the formulae from subsection IA.6.31 furnish 
the following expression for the complex Love numbers of a Kelvin- Voigt body: 

Ux) = ^77^ ^ ^ A ' (209) 

2(/ - 1) 1 + Ai {1 + Ty ix) 



Churkin (1998) employed the notation fc;(T) for what we call ki [t) . Our notations are more convenient in 
that they amplify the close analogy between the Love functions and the compliance function. 



51 



It then can be demonstrated, with aid of ( 16T|) . that the time-derivative of the corresponding Love 
function is 



2(/ 



(210) 



5(r) 



2(/-l) 1 + 

while the Love function itself has the form of 

3 1 



for Ty = 



2(/-l) A^C.r, 
3 1 



V 



1 - exp( - rC, 



e(r) 



(211) 



2(/ - 1) 1 + A 



l-exp(-rC) e(r 



where 



(212) 



Formulae (12101) may look confusing, in that exp( —tC,^) simply vanishes in the elastic limit, i.e., 
when — and — )■ oo . We however should not be misled by this mathematical artefact 
stemming from the nonanaliticity of the exponent function. Instead, we should keep in mind that 
a physical meaning is attributed not to the Love functions or their derivatives but to the results 
of the Love operator's action on realistic disturbances. For example, a Heaviside step potential 



Wi{R, r \ t') = WQ{t') 



(213) 



applied to a homogeneous Kelvin- Voigt spherical body will furnish, through relation ( l60l) . the 
following response of the potential: 



Uiir,t) = {-] / ki{t-t')We{t')dt' = (-) / ki{t-t')Wdt' 

Jt'=-oo \ / Jt'=0 



(214) 



w 



R 



l+l rT = t 



T J Jt = 

In the elastic limit, this becomes: 



ki (r) dr 



2{l - 1) 



1 - exp( -tC; 
1 + A, 



Uiir, t) ^ 



1 



2(/- 1) 1 + Ai \r 



R 



R 

r 



W 



W 



(215) 



O 1 

which reproduces the case described by the static Love number k^ = ^ 



An alternative way of getting fl215p would be to employ formulae (12 lip and (I59ap . 



C.2 A homogeneous Maxwell spherical body 



Using the formulae presented in the Appendix IA.6.41 and relying upon the Correspondence 
Principle, we write down the complex Love numbers for a Maxwell material as 



2(/-i: 



1 

1 + ^Af ^ X 



3 1 

2(/ - 1) 1 + A 



1 + 



1 + {1 + A,) T^ZX 



(216) 
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which corresponds, via ( pTl) . to 



k, r 



3 S{t) + 7; exp( -r-f ) 0(r) 



2{l-l] 



1 + A, 



(217) 



and 



1 + A, 



1 — exp( — 



1 + A, 



e(r) 



where 



7j 



(1 + A) 

A Heaviside step potential 

Wi{R, r*, t') = WQ{t') 
will, according to formula (l60|l . render the following response: 



(218) 



(219) 



(220) 



Ui{r, t) 



/+1 rt'=t 



kl {t-t')W Q{t') dt' 



t'= -oo 



l+l rt'=t 



kl {t-t')W dt' 



t'=0 



(221) 



l+l rT = t 



1 + A, 



1 - exp( -t7^ 



1 + A, 



R 



- ) ^ • 



In the elastic limit, we obtain: 
^ 7,->0 



t/^r, t) ^ 



2(/ - 1) 1 + A, 



W 



(222) 



which corresponds to the situation described by the static Love number k. 



2{i-i) 



C.3 A homogeneous Hohenemser-Prager (SAS) spherical body 

The Correspondence Principle, along with the formulae from subsection IA.6.5[ yields the following 
expression for the complex Love numbers of a Hohenemser-Prager (SAS) spherical body: 



kix) 

Combined with f l6T|) . this entails: 

■ , , 3 1 

k, (r) 



2(/-l) 1 + ^Xr, 

^ + 1 + ixr 



(223) 



y_ 

M 



2(/-l) 1 + A, ^ 

M 



A T — T 

5{t) + ' ^ 



1 + A, 



exp 



(224) 



and 



k,{r) 



1 _ Al. V ^A/ 



exp 



1 + A, 



T 



2(/-l) 



1 + ^/ 

M 



e(r) (225) 



53 



A Heaviside step potential 

Wi{R, r*, t') = WQ{t') (226) 
applied to a SAS spherical body will then result in the following variation of its potential: 

Uiir,t) = 

(j) kiit~t')Weit')dt'=i^\ k^^t-t')Wdt' = wij\ J kiir)dr 



2(1 - 1) 



\ V / 



1 — exp 



1 + A, 



- W 0(t) . (227) 
r / 



Within this model, the elastic limit is achieved by setting r^^ = , whence we obtain the 

O 1 

case described by the static Love number = ^ 1 + A ' Interestingly, the elastic regime is 

achieved even when these times are not zero. Their being equal to one another turns out to be 
sufficient. 

Repeating the above calculation for tidal disturbance WQ{—t') , we shall see that, after the 
tidal perturbation is removed, a tidally prestressed sphere regains its shape, the stress relaxing 

/ 1 + A 
at a rate proportional to exp I — j — -r-^ — t 

D The correspondence principle 
(elastic-viscoelastic analogy) 

D.l The correspondence principle, for nonrotating bodies 

While the static Love numbers depend on the static rigidity /i through ([3]), it is not immediately 
clear if a similar formula interconnects also ki{x) with /i(x) • To understand why and when 
the relation should hold, recall that formulae ([3]) originate from the solution of a boundary- value 
problem for a system incorporating two equations: 

a. = 2/iu. , (228a) 



<9(T^^ dp ^ d{W + U) 



the latter being simply the equation of equilibrium written for a static viscoelastic medium, 
in neglect of compressibility and heat conductivity. The notations a^^^ and u^^ stand for the 

deviatoric stress and strain, p = — ^ Sp S is the pressure (set to be nil in incompressible media), 
while W and U are the perturbing and perturbed potentials. By solving the system, one arrives 
at the static relation Ui = kiWi , with the customary static Love numbers ki expressed via p , 
R , and /i by ([3]) . 

Now let us write equation like f l228al - l228bl) for time-dependent deformation of a nonrotating 
body: 

§ = 2 /i U , (229a) 
p ii = V§ - Vp - V{W + U) (229b) 
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or, in terms of components: 

V = 2 /i , (230a) 

. ^ ^ _ _ djW + U) 

In the frequency domain, this will look: 

where a bar denotes a spectral component for all functions except fi - recall that /i is a spectral 
component not of the kernel /i(r) but of its time-derivative /i(r) . 

Unless the frequencies are extremely high, we can neglect the body-fixed acceleration term 
Uf^^ (x) in the second equation, in which case our system of equations for the spectral com- 
ponents will mimic fl228p . Thus we arrive at the so-called correspondence principle (also known 
as the elastic-vis coelastic analogy), which maps a solution of a linear viscoelastic boundary- value 
problem to a solution of a corresponding elastic problem with the same initial and boundary 
conditions. As a result, the algebraic equations for the Fourier (or Laplace) components of the 
strain and stress in the viscoelastic case mimic the equations connecting the strain and stress in 
the appropriate elastic problem. So the viscoelastic operational moduli p,{x) or J(x) obey the 
same algebraic relations as the elastic parameters /x or J . 

In the literature, there is no consensus on the authorship of this principle. For example, 
Haddad (1995) mistakenly attributes it to several authors who published in the 1950s and 1960s. 
In reality, the principle was pioneered almost a century earlier by Darwin (1879), for isotropic 
incompressible media. The principle was extended to more general types of media by Biot (1954, 
1958), who also pointed out some limitations of this principle. 



D.2 The correspondence principle, for rotating bodies 

Consider a body of mass Mprim , which is spinning at a rate uj and is also performing some 
orbital motion (for example, is orbiting, with its partner of mass Msec , around their mutual 
centre of mass). Relative to some inertial coordinate system, the centre of mass of the body is 
located at x^j^ , while a small parcel of its material is positioned at x . Relative to the centre of 
mass of the body, the parcel is located at r = x — x^j^j . The body being deformable, we can 
decompose r into its average value, Tq , and an instantaneous displacement u : 



X = x^j^j + tq + u . (232) 



Denote with D/Dt the time-derivative in the inertial frame. The symbol d/dt and its synonym, 
overdot, will be reserved for the time-derivative in the body frame, so dr^/dt = . Then 

Dr dr _ _ , D^r d'^f _ dr _ , _ _ , , 

= — + (jj X r and — - = — - + 2 lj x — + ujx{ujxr) + ujxr . 233) 

Dt dt Dt^ dt^ dt ^ ' ^ ' 

Together, the above formulae result in 

D'^x D'^x„., D'^r D'^x„., d'^r ^ ^ dr ^ _ , 

- ' - '^'''+ — + 2ivx — + cjx{i:uxr)+ujxf 



Df^ Dt^ Df^ Dt^ df^ dt 



+ ^ + 2c3x^ + c3x{c3xr)+dxr. (234) 
Ut^ dt'^ dt 
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The equation of motion for a small parcel of the body's material will read as 

P -j^ = V§ - Vp + F,,if + Fext , (235) 

— » 

where F^xt is the exterior gravity force per unit volume^ while Fgeif is the "interior" gravity 
force per unit volume, i.e., the self-force wherewith the rest of the body is acting upon the selected 
parcel of medium. Insertion of f l234p in f l235p furnishes: 

V§ - Vp + F.eif + F,xt . (236) 



P 



CM \ — '■lO— *■ — *■!— * /— * — *\ I — *■ — * 

+ n+ 2ci;x n+cjx [uj x r) + uj x r 



Dt^ 

At the same time, for the primary body as a whole, we can write: 

M,„m^^^ = J^^<^-td^r , (237) 

— » 

the integration being carried out over the volume V of the primary. (Recall that Fgxt is a force 
per unit volume.) Combined together, the above two equations will result in 



u + 2u3x u +UJ X {uj X r) + uj X r 



V§ - Vp + F.eif + Fext - -jf / Fext d^r . (238) 



For a spherically-symmetrical (not necessarily radially-homogeneous) body, the integral on the 
right-hand side clearly removes the Newtonian part of the force, leaving the harmonics intact: 

„ oo 

Fext - / Fext d'r = p J^^^i , (239) 

where the harmonics are given by 

W,{r,n = - {^) Pii^os^) , (240) 

r* being the vector pointing from the centre of mass of the primary to that of the secondary, 
and 7 being the angular separation between r and r * , subtended at the centre of mass of the 
primary. 

In reality, a tiny extra force F , the tidal force per unit volume, is left over due to the body 
being slightly distorted: 

» oo 

Fext - / Fext d^r = p + J" ■ (241) 

Here J-" is the density multiplied by the average tidal acceleration experienced by the body as a 
whole. In neglect of J-" , we arrive at 



P 



u + 2ujx u +UJ X {uj X r) + uj X r 



1=2 

Here, to each disturbing term of the exterior potential, Wi , corresponds a term Ui of the self- 
potential, the self- force thus being expanded into Fgei f = — S«^2 ■ 

Equation fl242p could as well have been derived in the body frame, where it would have 
assumed the same form. 

Denoting the tidal frequency with x ? we see that the terms on the left-hand side have the 
order of px^u , puxu , pu'^r , and pxur , correspondingly. In realistic situations, the first 
two terms, thus, can be neglected, and we end up with 

oo 

= V§ - Vp - p X V{Ui + Wi) - pGjx{ujxr) - pdxr , (243) 

1=2 

the term — Vp vanishing in an incompressible media. 
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D.3 The centripetal term and the zero-degree Love number 



The centripetal term in f l243p can be split into a purely radial part and a part that can be 
incorporated into the W2 term of the tide-raising potential, as was suggested by Love (1909, 

— > — *■ 

1911). Introducing the colatitude 0' through cos0' = ■ , we can write down the 
evident equality 



a; X (a; X r) = uj {uj ■ r) 



r uj 



1 / ^ ^2 
- (uj ■ r] 

2 ^ ^ 



1-2-2 

— u r 

2 



1—2-2 

-uj r 

2 



cos 



1 1 

The definition P2(cos0') = ^(Scos^c/)' — 1) easily renders: cos^0' = P2(cos0') + ^, 
whence: 



cD X (a; X 7^ = V 



- a; ^ r ^ [P2 (cos </> ') — 1 



We see that the centripetal force splits into a second-harmonic and purely-radial parts: 



p oj X {uj X r) 



^ a;VP2(cos, 



+ V 



P -2-2 

— u r 
3 



(244) 



(245) 



where we assume the body to be homogeneous. The second-harmonic part can be incorporated 
into the external potential. The response to this part will be proportional to the degree-2 Love 
number k2 . 

The purely radial part of the centripetal potential generates a radial deformation. This 
part of the potential is often ignored, the associated deformation being tacitly included into the 
equilibrium shape of the body. Compared to the main terms of the equation of motion, this 
radial term is of the order of 10~^ for the Earth, and is smaller for most other bodies. As the 
rotation variations of the Earth are of the order of 10~^ , this term leads to a tiny change in the 
geopotential and to an associated displacement of the order of a micrometer. 

However, for other rotators the situation may be different. For example, in Phobos, whose 
libration magnitude is large (about 1 degree), the radial term may cause an equipotential-surface 
variation of about 10 cm. This magnitude is large enough to be observed by future missions 
and should be studied in more detail^ The emergence of the purely radial deformation gives 
birth to the zero-degree Love number (Dahlen 1976, Matsuyama & Bills 2010). Using Dahlen's 
results, Yoder (1982, eqns 21 - 22) demonstrated that the contribution of the radial part of the 
centripetal potential to the change in mean motion of Phobos is about 3%, which is smaller 
than the uncertainty in our knowledge of Phobos' k2/Q . It should be mentioned, however, that 
the calculations by Dahlen (1976) and Matsuyama & Bills (2010) were performed for steady (or 
slowly changing) rotation, and not for libration. This means that Yoder's application of Dahlen's 
result to Phobos requires extra justification. 

What is important for us here is that the radial term does not interfere with the calculation 
of the Love number. Being independent of the longitude, this term generates no tidal torque 
either, provided the obliquity is neglected. 



D.4 The toroidal term 

The inertial term — p a; x r in the equation of motion 

— p l3 X r = p f X V(a; ■ r 



can be cast into the form 



(246) 



Tim Van Hoolst, private communication. 
^'^ Tim Van Hoolst, private communication. 
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whence we see that this term is of a toroidal type. Being almost nil for a despinning primary, 
this force becomes important for a librating object. 

In spherically-symmetric bodies, the toroidal force f l246p generates toroidal deformation only. 
This deformation produces neither radial uplifts nor variations of the gravitational potential. 
Hence its presence does not influence the expressions for the Love numbers associated with 
vertical displacement ( ) or the potential {k^). As this deformation yields no change in the 
gravitational potential of the tidally-perturbed body, there is no tidal torque associated with 
this deformation. Being divergence- free, this deformation entails no contraction or expansion 
either, i.e., it is purely shear. Still, this deformation contributes to dissipation. Besides, since 
the toroidal forcing results in the toroidal deformation, it can, in principle, be associated with a 
"toroidal" Love number. 

To estimate the dissipation caused by the toroidal rotational force, Yoder (1982) introduced 
an equivalent effective torque. He pointed out that this force becomes important when the 
magnitude of the physical libration is comparable to that of the optical libration. According to 
Ibid., the toroidal force contributes to the change of the mean motion of Phobos about 1.6%, 
which is less than the input from the purely radial part. 



E The Andrade and Maxwell models at different frequen- 
cies 

E.l Response of a sample obeying the Andrade model 

Within the Andrade model, the tangent of the phase lag demonstrates the so-called "elbow 
dependence". At high frequencies, the tangent of the lag obeys a power law with an exponent 
equal to —a , where < a < 1 . At low frequencies, the tangent of the lag once again obeys a 
power law, this time though with an exponent — (1 — a) . This model fits well the behaviour 
of ices, metals, silicate rocks, and many other materials. 

However the applicability of the Andrade law may depend upon the intensity of the load 
and, accordingly, upon the damping mechanisms involved. Situations are known, when, at low 
frequencies, anelasticity becomes much less efficient than viscosity. In these cases, the Andrade 
model approaches, at low frequencies, the Maxwell model. 



E.1.1 The high-frequency band 

At high frequencies, expression fl89b|) gets simplified. In the numerator, the term with 



1 



dominates: z ^ z ^C, which is equivalent to z ^ " . In the denominator, the constant 

term dominates: 1 ^ , or simply: z ^ 1 . To know which of the two conditions, 

1 

z ^ " or z 1 , is stronger, we recall that at high frequencies anelasticity beats viscosity. 
So the a— term in (185|) is large enough. In other words, the Andrade timescale should be 
smaller (or, at least, not much higher) than the viscoelastic time r^^ . Accordingly, at high 

frequencies, ( is smaller (or, at least, not much higher) than unity. Hence, within the high- 

1 

frequency band, either the condition ^ ^ 1 is stronger than z ^ " or the two conditions 
are about equivalent. This, along with (!90|) and (I9T]) enables us to write: 

tan5 ^ (xrj sin (— j r(a + 1) for x » r'' = t^'C' ■ (247) 

The tangent being small, the expression for sin 6 looks identical: 

sin5 ^ ixrS" (t^) ^(^ + 1) for X » r;' = r-JC' ■ (248) 
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E.1.2 The intermediate region 



In the intermediate region, the behaviour of the phase lag 6 depends upon the frequency- 
dependence of ( . For example, if there happens to exist an interval of frequencies over which 

1 

the conditions 1 ^ 2; ^ " ^-^e obeyed, then over this interval we shall have: 1 ^ 
and z'" ^ C^~^ • Applying these inequalities to fl89bl) . we see that over such an interval of 



frequencies tan 5 will behave as z tan ^'^^ 



E.1.3 The low-frequency band 

At low frequencies, the term z^^ ( becomes leading in the numerator of ([89b]): z'"' ^ C , 

1 

which requires z ^ " . In the denominator, the term with z'" becomes the largest: 1 <^ 
z~°' , whence z <^ 1 . Since at low frequencies the viscous term in ( 185|) is larger than the 
anelastic term, we expect that for these frequencies ( is larger (at least, not much smaller) than 

unity. Thence the condition z ^ 1 becomes sufficient. Its fulfilment ensures the fulfilment of 

1 

z <^ " . Thus we state: 

tan5 ^ ixrS'^''^^ T^^^f- " for x « r'' = r;^' C • (249) 

cosf^j r(a + l) 

The appropriate expression for sin 6 will be: 

sin5 ^ 1 - 0{{xr,y^'-''H-') for X « r;' = r;^ , (250) 



It would be important to emphasise that the threshold r^^ = r^^^ C,^^ standing in fl247p and 
( I248p is different from the threshold = showing up in f l249p and fl250p . even though 

these two thresholds are given by the same expression. The reason for this is that the timescales 

and r^^ are not fixed constants. While the Maxwell time is likely to be a very slow function 
of the frequency, the Andrade time may undergo a faster change over the transitional region: 
must be larger than r^^^ at low frequencies (so anelasticity yields to viscosity), and must become 
shorter than or of the order of r^,^ at high frequencies (so anelasticity becomes stronger). This 
way, the threshold standing in 02491 - [2^0]) is lower than the threshold standing in ( 12471 
- I248P . The gap between these thresholds is the region intermediate between the two pronounced 
power laws ([217]) and fl2l9D . 

E.1.4 The low-frequency band: a special case, the Maxwell model 

Suppose that, below some threshold Xo ; anelasticity quickly becomes much less efficient than 
viscosity. This would imply a steep increase of ( (equivalently, of ) at low frequencies. Then, 
in (]89bp . we shall have: 1 ^ z~°' and 2;"" ^ C^~^ ■ This means that, for frequencies below Xq > 
the tangent will behave as 

tan5 ^ = (xv)~^ for X < Xo • (251) 

the well-known viscous scaling law for the lag. 

The study of ices and minerals under weak loads (Castillo-Rogez et al. 2011, Castillo-Rogez 
2011) has not shown such an abrupt vanishing of anelasticity. However, Karato & Spetzler (1990) 
point out that this should be happening in the Earth's mantle, where the loads are much higher 
and anelasticity is caused by unpinning of dislocations. 
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E.2 The behaviour of sine;(x) = — Tm [/^/(x)] 

within the Andrade and Maxwell models 

As we explained in subsection 4.1, products ^^(X/mpg) ^^^^^(^impij) enter the Impq term of the 
Darwin-Kaula series for the tidal potential, force, and torque. Hence the importance to know 
the behaviour of these products as functions of the tidal frequency X/mpg • 



E.2.1 Prefatory algebra 

It ensues from ([2HD that 

3 (7^e[J(x)])' + (Xm[J(x)])' \ A, J 'Re{J{x)\ +iA^JXm[J{x 



kiix) 



2 {I 



whence 



(7^e[J(x)] + A,jy + (Xm[J(x)]) 
_ 3 - Ai Jim [J{x 



(252) 



2 (/ - 1) ( 7^e [Jix)] + A,jY + {Im [J(x)] )' 



(253) 



J = J(0) = = l/yu(0) being the unrelaxed compliance (the inverse of the unrelaxed shear 
modulus yu). For an Andrade material, the compliance J in the frequency domain is rendered 
by (186|1 . Its imaginary and real parts are given by ([87|-[88]). It is then easier to rewrite (I253P as 

\kix)\ sine^(x) = 



3 A, 



C + sin (^^) T{a + 1) 



2(/-l) 
where 

and 



Ai + 1 + z-^ cos ( ^ j T{a + 1)J + z-^ + -2"" sin ) 1(1 + a) 



2 ' 



(254) 



(255) 



(256) 



For /3 — , i.e., for — )■ oo , f |254p coincides with the appropriate expression for a spherical 
Maxwell body. 

E.2. 2 The high-frequency band 

Within the upper band, the term with z'"' dominates the numerator, while A^ dominates the 

1 

denominator. The domination of z~°' in the numerator requires that z ':$> " , which is the 



same as x ^ ^ C ^ ° • The domination of A, in the denominator requires: z ^ A 



-l/a 



which 



is the same asx^''' C A 



l/a 



It also demands that ( z ^ <^ A^ , which is: x ^ t^,/ A^ ^ . 

■ l/a 



For realistic values of A^ (say, 10^) and a (say, 0.25), we have: A ~ 10 . At high 
frequencies, anelasticity beats viscosity, so ( is less than unity (or, at least, is not much larger 



than unity). On these grounds, the requirement x ^ '''m^ C 

.-1,-1 



1-a 



fulfilment guarantees that of both x^ '^m 

3 



is the strongest here. Its 
and X ^ ^A/M/ • Thus we have: 



\ki{x)\ sine;(x) 



2(/-l) A 



1 /an 
sm 



r(a + l) (r,,x; 



for 



_i ^ 1— a 

X » rJC 



(257a) 
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E.2.3 The intermediate band 

Within the intermediate band, the term ( takes over in the numerator, while still dom- 

1 

inates in the denominator. The domination of C, in the numerator implies that z <^C^^ " , 

a 

which is equivalent to x ^ C ^ " • The domination of in the denominator requires 

— l/a 

X ^ r~ C A^ and x ^ '^m A" ' j^^^ above. 

As we are considering the band where viscosity takes over anelasticity, we may expect that 

here ( is about or, likely, larger than unity. Taken the large value of A^ , we see that the condition 

11 11^"'"/'^ 
Y ^ A~ is stronger. Its fulfilment guarantees the fulfilment of v ^ C A, . This 

way, we obtain: 

\kix)\ sme,(x) ^ ^^y^ (r,,x)-' for r"^ C ^ » X » r^^A"' • (257b) 

E.2.4 The low-frequency band 

For frequencies lower than t^^^A^^ the Andrade model renders the same frequency-dependency 
as that given (at frequencies below ) by the Maxwell model: 

sine,(x) ^ 2{l^-l) ^ A"' > ^ • (257c) 

E.2.5 Interpretation 

Formulae fl257al) and fl257bl) render a frequency-dependence mimicking that of | J (x) | sin 5{x) 
= — Xm [ J(x)] in the high- and low-frequency bands. This can be seen from comparing f l257ap 
and (l257bl) with (12481) . 

In contrast, f l257cp reveals a peculiar feature inherent in the tidal lagging, and absent in the 
lagging in a sample. 

For terrestrial bodies, the condition Aj^ ^ x puts the values of x below 10~^°i/2, give 
or take several orders of magnitude. Hence |A;;(x)| sinej(x) follows the linear scaling law fl257cp 
only in an extremely close vicinity of the commensurability where the frequency x vanishes. 
Nonetheless it is very important that |A;2(x)| sine(x) first reaches a finite maximum and then 
decreases continuously and vanishes, as the frequency goes to zero. This confirms that neither 
the tidal torque nor the tidal force becomes infinite in resonances. 



F The behaviour of ki{x) = |^/(x)l ^he 
hmit of vanishing tidal frequency x ? 
within the Andrade and Maxwell models 



From (168!) . we obtain: 

\h{x)\' ' ^ 



Jix)\ 



2{l-l)J \A,J+ J{x)\ V2(/-2) 



7^e 






J ix)\ 












J ix)_ 


r 



Bringing in expressions for the imaginary and real parts of the compliance, and introducing 
notations 

E = 7]-^ , 5 = /3 sin (^^) r(l + a) , D = /3 cos ( ^ ) 1(1 + a) 



(259) 



(258) 
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we can write: 

\h{x)? = 

and 

\h{x)\-' = 

whence 

\Hx)\ - 

and 

\Hx)\-' 





l)J 


f2{l - 


l)\ 


V 3 


) 


3 




2{l- 


1) 


2{l - 


-1) 



1 ^ x'"" ' ' — x' + oix'-'n 



E^ 



D 



E^ 



the expansions being vahd for xJ/E <^ 1 + A ; i-^-, for X'^m ^ ^ + A 
Rewriting fl68|) as 



(260) 



,(261) 



(262) 



(263) 



Hx) 



Tie 



Hx) 



i Xm 



Hx) 



Tie 



Hx) 



A^ J — i Xm 



Hx) 



2(/-l) 
we extract its real part: 



Tie 



Hx) 



+ A + (Xm 



Hx) 



(264) 



7^e 



Hx) 



Xm 



Hx) 



A^JTZe J ix) 



Tie 



Hx) 



+ A -^1 + {Irn 



Hx) 



2(1-1) 



Tie 



Hx) 



+ A, J 



Tie 



Hx) 



Insertion of the expressions for IZe J (x) 



+ A,J 
and Xm 



Xm 



Hx) 



(265) 



Hx) 



ne [hix)] 



2(/-l) 



D 

l-A,J-x 



A J 



into the latter formula entails: 

. (266) 



Expressions (12631) and (12661) enable us to write down the cosine of the shape lag: 



cos ei 



Ux)\ 



E 



x' + Oix 



2 A(^mX)' + 0{x'-'n -(267) 



Comparing this expression with (I262p . we see that, for the Andrade {a ^ 0) model, the 
evolution of ki[x) = |^2(x)l the hmit of small x > unfortunately, cannot be approximated with 
a convenient expression ki{x) ~ ki{0) cose(x) , which is valid for simpler models (like the one 
of Kelvin- Voigt or SAS). 

However, for the Maxwell model (/3 = 0) expression (12621) becomes: 



Ux)\ 



2(/- 1) [ 



(J + A,J/2)x' + 0ix') 



2(1-1) 



[l-A,{l + Aj2){r,,x)' + 0{x')] m 
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which can be written as 



\h{x)\ 



2{l-l) 



1 

2 



E 



2{l-l] 



\A]{r,,xf + 0{x')\ ,(269) 



insofar as ^ 1 . Comparing this with f l267p . we see that, for small terrestrial moons and 
planets (but not for superearths whose Ai is small), the following convenient approximation is 
valid, provided the Maxwell model is employed: 



h{x) ~ ^/(O) cose(x) , for x'T'm ^ ^ + A where > 1 



(270) 



G The eccentricity functions 



In our development, we take into account the expansions for G^^^^ 
tricity, keeping the terms up to 



e) over the powers of eccen- 
inclusively. The table of the eccentricity functions presented 



in the book by Kaula (1966) is not sufficient for our purposes, because some of the Gi^^^e) 
functions in that table are given with lower precision. For example, the term is missing in the 
approximation for G^qq (e) . Besides, that table omits several functions which are of order . 
So here we provide a more comprehensive table. The table is based on the information borrowed 
from Cayley (1861) who tabulated various expansions employed in astronomy. Among those, 
were series 





[cos]* 


cos 




[sin]* 


sin 



E 



cos 
sin 



iM 



(271) 



V and M. signifying the true and mean anomalies, while [cos]* and [sin]* denoting the co- 
efficients tabulated by Cayley. These coefficients are polynomials in the eccentricity. Cayley's 
integer indices i , j are connected with Kaula' s integers I , p , q via 



I -2p = J 



I - 2p + q 



(272) 



With the latter equalities kept in mind, the eccentricity functions, for i > , are related to 
Cayley's coefficients by 



cos ■+ sin 



for i > . 



(273) 



To obtain the eccentricity functions for i < , one has to keep in mind that [ cos 
while [sin]* = — [sin]*. It is then possible to demonstrate that 



cos 



cos — sm 



for i < 



(274) 
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Then the following expressions, for 1 = 2, ensue from Cayley's tables: 

= ^20 -10 (e) = G^o -9(e) = G^o -8(e) = , 



G20 -11 (e 
G20 -ri^ 
G20 -ei^ 
G20 -5(6 
G 
G 
G 

G20 -1' 



•20 -4^ 



'20 -3* 



■^20 -2^ 



15625 7 
e 



129024 
4 



45 



1280 
1 



81 
2048 



7 



— + e" 

24 240 



1 o 11 , 313 7 

— e + e + e 

48 768 30720 



= , 



e + — e-" - 

2 16 384 



e' - 



143 
18432 



(275a) 
(275b) 
(275c) 
(275d) 
(275e) 
(275f) 
(275g) 
(275h) 



G'2oo(e) 
<^2oi(e) 
<^202(e) 



1 e^H 

2 16 



7 123 o 

- e - e"^ + 

2 16 128 



288 

489 . 1763 



2048 



17 2 115 4 601 

— e - e + 

2 6 48 



(275i) 
(275j) 
(275k) 



^ , , 845 o 32525 , 208225 
2°=^^ ^ 48 768 6144 



(2751) 



G'204(e) 
G'205(e) 

G'206(e) 

G'207(e) 



533 . 13827 



16 

228347 
3840 

73369 6 
720 ^ 

12144273 
71680 



160 



, 3071075 7 
- — — 



18432 



(275m) 
(275n) 
(275o) 
(275p) 



the other values of q generating polynomials G^q^ (e) whose leading terms are of order and 
higher. 

Since in our study we intend to employ the squares of these functions, with terms up to 
only, then we may completely omit the eccentricity functions with |g| > 4. In our approximation. 
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the squares of the eccentricity functions will look: 



(e) = + O(e^) , (276a) 

20 -3 ^ ' 2304 ' V ; ' \ / 

Gl_,{e) = , (276b) 

(e) = ^ - 1 + + 0(e«) , (276c) 

20 -1 ^ ^ 4 16 768 y > ^ y > 

GMe) = 1 - 5e^ + - - — + 0(e«) , (276d) 
200 ^ ^ 8 36 

= ^e^-if5,«+ O(e^) , (276f) 

the squares of the others being of the order of or higher. 

Be mindful that, for / = 2 we considered only the functions with p = . This is dictated 
by the fact that the inclination functions F^^^ = F^^^ are of order i (and, accordingly, their 
squares and cross-products are of order ) for all the values of p except zero. 

For / = 3 , the situation changes. The inclination functions Fi^p = -F310 , -^312 , -^313 , -^320 ) 
-^321 , -^322 , -P323 , -^331 , -^332 , -^333 are of the order 0{i) or higher. The terms containing the 
squares or cross-products of the these functions may thus be omitted. (Specifically, by neglecting 
the cross-terms we get rid of the mixed-period part of the / = 3 component.) What is left is the 
terms with Imp = 311 and Imp — 330. These terms contain the squares of functions 

F3ii(«) - -\ + 0{e) and F33o(^) = 15 + Oil") , (277) 

accordingly. Prom here, we see that, for Z = 3 , we shall need to employ the eccentricity functions 
Gi,,{e) = G,,^{e) and G,^^{e) = G,,^{e) . 

The following expressions, for / = 3 and p = , ensue from Cayley's tables: 

^30-11(6) = G,,_,,{e) = G,,_,{e) = G,, _,{e) = , (278a) 

G,,.,{e) = . (278b) 

G3o-6(e) = ' (278c) 

1 1 3 

G30 -.(e) = Bo e= + liio • P78d) 
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^30 -3(^) 


= , 


(278f) 


^30 -2(^) 


= ^ + ^ + , 
8 48 3072 


(278g) 




= - e + - e-^ - — e' , 

4 48 288 


(278h) 


^300 (^) 


. 2 ^ 423 4 125 6 
= ^ + 64 ^ 64 ^ ' 


(278i) 


G'30l(e) 


= 5e-22e3+60^^-9^^ 
24 9 


(278j) 


^302 (^) 


127 9 3065 . 243805 « 

= e - H , 

8 48 3072 


(278k) 


^303 (^) 


163 q 2577 r; 1089 7 

= e"^ - H , 

4 16 5 ' 


(2781) 


^304(6) 


35413 4 709471 g 
~ 384 ^ 1920 ^ ' 


(278m) 


^305(6) 


23029 . 35614 7 
120 45 


(278n) 



385095 
1024 

44377 
-63-^ 



(278o 
(278p 



the other values of q generating polynomials G^^^ (e) and G^^^ (e) , whose leading terms are of 
order and higher. 

The squares of some these functions, will read, up to terms inclusively, as: 



(e) 

30-3^ ^ 



(e) 

30-2^ ^ 



G' (e) 

30 -1 ^ ^ 



G' (e) 

300 ^ ^ 



G' (e) 

301 ^ ^ 



G' (e) 

302 ^ ^ 







64 192 ^ 



1 - 12 + 
25 - 220 + 



1575 4 2663 ^ ^,0, 



32 32 

8843 



+ 0(e«) , 



16129 4 389255 « 
e — 



Gl(e) = 



303 



64 

26569 
16 



192 



12 

+ 0(e«) , 



e' + 0{e') , 



(279a 
(279b 
(279c 
(279d 
(279e 
(279f 
(279g: 



the squares of the others being of the order of or higher. 
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Finally, we write down the expressions for the eccentricity functions with / = 3 and p = 1: 











16337 7 

— c 
2240 


('280h'l 


'-'si -&\^) 


48203 6 
~ 9216 ^ ' 






899 , 2441 7 
— c 1 c 
240 480 


(2R0r\) 


G (p'] 

'-^31 -A\^) 


343 4 2819 

— r 1 r" 

128 640 


f280p'l 


G (p) 


23 89 , 5663 . 

— r \ r \ r 

12 24 960 


f280f'l 


<^31 -2(e) 


11 2 49 4 15665 
lb oO/z 


(280g) 


G^i _i(e) 


5 35 , 105 7 

= e H e'^ H H , 

2 8 16 ' 


(280h) 


G'siole) 


239 /I 3323 ^ 

= 1 + 2 , 

64 576 


(280i) 


G3n(e) 


11 245 , 463 7 
= 3e+ — e^+ — 6^+ — 
4 48 64 


(280j) 


G Cp^ 


53 2 39 4 7041 . 

— r \ r \ r 

8 16 1024 


I ZiOUJv J 


G (p^ 

^313V^^ 


77 , 25 , 4751 , 
6 48 480 


('2801'! 


G (p^ 

^314 V^^; 


2955 4 3463 ^ 
128 " 384 


('OSOrri'l 


G (p^ 


3167 , 8999 , 
80 " 320 


/'OSOn'l 




3024637 6 
46080 ^ ' 


(280o) 




178331 7 
1680 ^ ' 


(280p) 



and the squares: 



G' (e) 

31-3^ ' 



G (e) 

31 -2 ^ ' 



G' (e) 

31 -1 ^ ' 



(e) 

310^ ' 



529 
144 



+ 5 + 15 + 0(6^) , 



1 + 4 + 



367 



7625 
^88" 



(281a) 
(281b) 
(281c) 
(281d) 



G' (e) 

311 ^ ^ 



( 

312 



G' (e) 

313 ^ ' 



9 + ^ + ^ + 0(e«) 
2 16 ^ ^ 



2809 
64 

5929 
^6~ 



2067 



+ 0(e« 



+ 0(e« 



the squares of the other functions from this set being of the order or higher. 



(281e) 
(281f) 
(281g) 



H The / = 2 and / = 3 terms of the secular part of the 
torque 

H.l The / = 2 terms of the secular torque 



Extracting the / = 2 input from flllSp . we recall that only the [Impq) = (220g) terms matter. 
Out of these terms, we need only the ones up to . These are the terms with |g| < 3 . They 
are given by formulae f l276p from Appendix [Gl Employing those formulae, we arrive at 



r 



1=2 



T 



(«mp)=(220) 



+ 0{^\) = -GMi^R'a-^ ^L^'^ + ^^"'^^ + (282a) 



g=-3 



220g 



-GMl^R'a- 



ko sine 

23 4 220 -3 



1 9 -^4 fi \ J 

- e — — e + e k2 sine„„„ , 

4 16 768 / 220 -1 
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155 
~36 



+ ( 1 - 5 + — - — eM fca sin + — 



"2200 



49 o 861 



4 



16 



21975 . 

e Ko sine 

256 / 2201 



+ 



289 



1955 



, . 714025 6 , . 
k-} sm e H e ko sm e 

2304 2203 



'2202 



+ 0(e«e) + 0(^2 e) 



(282b) 



where the absolute error 0(e*e) has emerged because of our neglect of terms with |g| > 4 , 
while the absolute error 0{i^ e) came into being after the truncation of terms with p > 1 . 
Recalling expression f ll09bp . we can rewrite (1282bp in a form indicating explicitly at which 



resonance each term changes its sign. To this end, each h sine will be rewritten as: 
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ki sin |e, I sm 

' Impq ' 



{I — 2p + q) n — mO 
1 



. This will render: 



1=2 2 



e ko sin I e I sgn ( — n — 2 

2304 I 220 -3 I 



1 1 13 

^ + e" 1 fco sin |e | sgn ( — 2 6' 

7g8 I ^ I 220 -1 ' ' 



+ I r 16 



It ^2 4 155 h \ 7 • I 

+ 1 — 5e + — e — e /co sm e | sgn ( n 



'49 2 861 4 21975 ^ , , , , 



289 4 1955 « , , . I I . 

e /C2 sm e„„„„ sen { 2n — 6 



+ ' — " 6 



2202 I 



714025 

+ — fco sinle I sgn ( 5n — 3 6* 

z I 2203 I ^ \ 



2304 



+ 0(e«e) + 0(z2e) , (283) 



H.2 The 1 = 3, m = 1 terms of the secular torque 

Getting the Z = 3 , m = 1 input from (11131) and leaving in it only the terms up to , we obtain, 
with aid of formulae ( 1276^ from Appendix O, the following expression: 



r = -GMlR'^a-^y^G^ (e)kssme +0(e^e) 

amp)=(311) 8 ^ 3ig^ ^ 31ig ^ ^ 

g=-3 



(284a) 



G a 



529 6 , . 121 4 539 y , , . 

e fccj sme + e + e fcq sme 

144 ^ -3 I 64 64 / ^ 311 -2 



+ (e^ + 5e^ + 15e^) k, sine^^^ -i + + ^ + ^ + ^ 



1 ^3 sine 



3110 



. 33 4 611 A , /2809 , 2067 . , , 

9 H + sm e + H A;, sm e 

2 16 / ^ 3111 ^ V 64 64 I ^ 3112 



5929 fj, 
+ A^s sme 

3d 



3113 



+ O(e^) 



(284b) 
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With the signs depicted exphcitly, this will look: 



(imp)=(311) 8 ^""^ 



529 
144 



ks sinle^^^ _^\ sgn (^-2n - 9^ 



12X ^ 539 (s \ 7 • I I / n 

+ ' "61 ^ i ''''''311 -2I ( - ^ 



+ (e^ + 5 + 15 e^) /cs sin \€^^^ _ J sgn (| - ^ ) 



367 , 7625 



+ 1 + 4e^+ — + 



32 288 



e^^ ks sinle^^^^l sgn - 



(33 611 ^ 

Qe" + y + — eM /cs sinjeg^^J sgn (^2n - 9^ 



+ 



/2809 4 , 2067 



^4 ^ 2067 sin|e3^^^| sgn [^n - 9^ 



+ 



5929 



e6A;3 sin 163^^3 I sgn (4n - ^) 



H.3 The / = 3 , m = 3 terms of the secular torque 

The second relevant group of terms with I — 3 will read: 
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15 

(/mp)=(330) 8 



q=-3 



GMl^R'a-' E ^3og^^) ^3 sine33^^ + 0(e«e 



3305 



64 192 



— + 7-^ 1 /cs sine 



330 -2 



+ I 5e^+ U3 sine_ +(l-12e" + 



2 48 



330 -1 



1575 , 2663 



32 32 



8843 



/ 16129 , 389255 



+ |26e'- 220 e'+^e»jfc3sm.33., + V 64 " 192 



26569 fj, 
+ — TT^e kz sine 
16 



3303 



70 



or, with the signs shown exphcitly: 
15 



r 



(/mp)=(330) 8 



+ e 



(e^- \e'+ ge^) k, s^^\e^,, _J sgn (-^) 



_ 2 1575 4 2663 



I sin |e I sgn [ n — 9] 

J I 3300 ' ° V / 



+ 25 - 220 + 



8843 
1^ 



e^^ /C3 sinle^g^J sgn ^2n - 6^^ 



/ 16129 4 389255 
+ I _ . e 



V 64 



192 



e') ^3 Sin 1633^2 I sgn (3n - ^) 



26569 6, • I I n 
+ e /cq sm e sgn [ An — 

15 I 3303 1 o \^ 



+ 0{e'e) 



(287) 



I The / = 2 and / = 3 terms of the short-period part of 
the torque 

The short-period part of the torque may be approximated with terms of degrees 2 and 3: 
f ^ %^, + ^3 + Oie{R/ar) 



T + 

(imp) =(220) 



r + r 

{lmp)=(3U) (imp)=(330) _ 



+ 0{et^) + 0{e{R/a)^) , (288) 



where 
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(/mp)=(220) 



3 3 



SGM^^^R a-^ Yl 5] ^^20,(6) G2o,(e) { cos[M{q-j)] A;^ sine^^o? 



?=-3 j=-3 



sin[A^(g-j)] ^2 cose^^Q^ } + 0(^2 e) + 0(e^ e) , (289a) 



-3GMli?'a-6 ^ J2G2o,{e)G2oj{e)k^sm[M{q-j)]+Oiz'e)+0{ee) , (289b) 



3 3 



-3 j=-3 

j < g 
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W(3n) = fl jl G,,,{e) G,,,{e) {cos [A^ (g - ,) ] h^u^e^^^^ 

9=-3 j=-3 

- siYi[M{q- j)] ^2 cose^„ | + 0(i^ e) + 0(e^ e) , (290a) 



o 3 3 

3 _ , 2 _5 



GM,,,i?a-6^ Y,G?M{e)G^i,{e)k^sin[M{q-3)] + 0{ee) + 0{ee) , (290b) 

g=-3 j=-3 



1 r 3 3 

Vp)=(330) = -^GM^^^R a-' J2 J2 G-soqie) Gsojie) cos[A^(g-j)] h ^^^^ssoq 

q=-3 j=-3 



- sm[M{q-j)] cose^g^^ }+ 0{f e) + O(e^e) , (291a) 



^^GMl^Ra-' J2 i2G30q{e)Gso,{e)k,sm[M{q-j)]+0{t'e) + 0{ee) , (291b) 



4 



j < q 



the expressions for the eccentricity functions being provided in Appendix [Gl The overall nu- 
merical factors in (12891 - I29ip are twice the numerical factors in (11141) , because in (12891 - I29ip 
we have j < q and not j ^ q- The right-hand sides of (12891 - 12911) contain 0{ee) in- 
stead of O(e'^e) , because at the final step we approximated cos[A^(g — j)] k^ sine^^^^ — 
sin[A^(g — j)] ki cose^^^^ simply with — sin[A^(g— j)] k^ . Doing so, we replaced the 
cosine with unity, because the entire Darwin-Kaula formalism is a linear approximation in the 
lags. We also neglected k^ sin e^^^^ and kept only the leading term with k^ . This neglect would 
be illegitimate in the secular part of the torque, but is probably acceptable in the purely short- 
period part, because the latter part has a zero average and therefore should be regarded as a 
small correction even in its leading order. The latter circumstance also will justify approximation 
of ki = kiix) with fc^(O) in (1289) - [29T]) . 
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